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OPTIMAL GUIDANCE LAW DEVELOPMENT 
FOR AN ADVANCED LAUNCH SYSTEM 


Anthony J. Calise* and Martin S. K. Leung** 

Georgia Institute of Technology, GA 30332 

SUMMARY 

The objective of this research effort was to develop a real-time guidance approach 
for launch vehicles ascent to orbit injection. Various analytical approaches combined 
with a variety of model order and model complexity reduction have been investigated. 
Singular perturbation methods were first attempted, and found to be unsatisfactory. The 
second approach based on regular perturbation analysis was subsequently investigated. It 
also fails because the aerodynamic effects (ignored in the zero order solution) are too 
large to be treated as perturbations. Therefore, the study demonstrates that perturbation 
methods alone (both regular and singular perturbations) are inadequate for use in 
developing a guidance algorithm for the atmospheric flight phase of a launch vehicle. 

During a second phase of the research effort, a hybrid analytic/numerical 
approach was developed and evaluated. The approach combines the numerical method of 
collocation and the analytical method of regular perturbations. The concept of choosing 
intelligent interpolating functions is also introduced. Regular perturbation analysis 
allows the use of a crude representation for the collocation solution, and intelligent 
interpolating functions further reduce the number of elements without sacrificing the 
approximation accuracy. As a result, the combined method forms a powerful tool for 
solving real-time optimal control problems. Details of the approach are illustrated in a 
fourth order nonlinear example. The hybrid approach is then applied to the launch 
vehicle problem. The collocation solution is derived from a bilinear tangent steering law, 
and results in a guidance solution for the entire flight regime, that includes both 
atmospheric and exoatmospheric flight phases. Assessment of performance and 
reliability are demonstrated through closed loop simulations. The hybrid guidance 
approach delivers over 99.9% of optimal performance and orbit injection accuracy while 
the control computation is completed in tenths of a second on a SPARCstation 1. Wind 
shear effects and a control constraint are also addressed. 


* Professor, School of Aerospace Engineering. 
** Graduate Research Assistant. 
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A second effort that paralleled this work under the same grant number was lead by 
Dr. Dewey Hodges, of the School of Aerospace Engineering at Georgia Tech. This work 
has been documented under a separate contractor report. 
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SECTION I 


INTRODUCTION 

The objective of the Advanced Launch System (ALS) program is to develop an 
unmanned, all-weather launch system for placing large payloads (100,0001b - 150,0001b) 
into a low Earth orbit at a fraction of present cost. Part of the guidance requirement is to 
realize an efficient algorithm for solving the launch vehicle ascent trajectory problem. 

1.1 Background 

To date, first stage guidance has been realized in open loop form. The vehicle is 
typically guided by using a pre-stored steering program. The steering program is calculated 
as a part of pre-launched preparation to account for structural loads from aerodynamic 
forces and from atmospheric disturbances such as wind shear. Typically it involves flying 
with nearly zero angle of attack, and performing a gravity turn [1]. Near zero angle of 
attack is employed to avoid creating excessive aerodynamic bending moments, which is 
proportional to the product of angle of attack and dynamic pressure. Guidance for the 
second stage and any subsequent stages is closed loop, employing various approaches. 
The Saturn V vehicle uses an Iterative Guidance Mode (IGM) [2], and the Space Shuttle 
employs Powered Explicit Guidance (PEG) [3]. These are retargeting schemes because 
the guidance commands are recalculated at each update cycle using the current vehicle's 
position and velocity vectors as the initial conditions for the optimization process. 
Traditional Guidance Solution Methods 

Traditional launch vehicle guidance may involve either two or three different phases 
[1 - 3]. The first is an open loop guidance phase for the atmospheric portion of flight 
which typically flies with a non-optimal piecewise linear attitude program. The second is a 
closed loop guidance phase for the exoatmospheric portion of flight. This has an analytic 
solution under certain assumptions. Then a third closed loop phase is possibly required 
when the vehicle is approaching orbital conditions for final precision orbit injection. 

Numerical approaches to optimal guidance typically employ either nonlinear 
programming [4 - 9] or multiple shooting [10]. In a direct method formulation such as 
nonlinear programming, the optimization problem is transformed into a parametric 
optimization problem. The unknown control profile is parameterized with undetermined 
coefficients of typically piecewise linear polynomials. The states are considered as 
functions of the control through the differential equations of dynamics. Constraints, if any, 
are enforced discretely along the trajectory, typically at a finite number of nodal points of 
the parameterized control. So the original infinite dimensional problem is approximated by 



a finite dimensional problem in the reduced space of the control parameters, and gradient 
techniques are used to search for a solution that optimizes the performance index. In [8], 
Hargraves and Paris have combined the nonlinear programming method with collocation by 
approximating all the state and control histories with piecewise smooth functions, thus 
avoiding any integration process. Similar to the collocation method, Pamadi [9] has used 
splines as function of velocity to approximate the altitude profile and applied an 
optimization algorithm to determine the unknown coefficients of the splines. To be useful 
as a feedback guidance solution, it is essential that these approaches converge quickly and 
reliably at each instant the solution is updated during the flight 

On the other hand, multiple shooting is a technique used in indirect methods. 
Instead of evaluating the performance index directly, optimization is achieved by satisfying 
a set of necessary conditions which are expressed in the form of a Two-Point Boundary 
Value Problem (TPBVP). For a constrained case, this may lead to a Multi-Point Boundary 
Value Problem (MPBVP), for which a guess of the switching structure is required. To 
reduce the sensitivity to an initial guess of the solution, piecewise integration or multiple 
shooting is used. Instead of integrating for the complete trajectory starting from one set of 
initial conditions, the trajectory is divided into intervals and integration is performed 
separately from different sets of initial conditions for each interval. Then the boundary 
conditions and continuity conditions (or jump conditions in the case of state constraints or 
discontinuous dynamics) between intervals are enforced. A relaxed Newton’s method [11] 
is typically used to iterate for a solution. Though the indirect method produces extremely 
accurate results, it involves complicated programming in formulating the costates 
differential equations and the control structure. The process is also complicated by the 
requirement to provide an initial guess for both costate and state variables. On the contrary, 
nonlinear programming is relatively simple to formulate. The method does not require the 
use of costate variables or a knowledge of switching structure. In practice, it is favored 
over indirect methods for solving optimization problems in general purpose programs. 

Due to the intensive computation requirements, direct and indirect methods are used 
only to generate off-line solutions for analysis purposes or to provide a first stage open 
loop guidance program. To compensate for using an open loop approach during the first 
stage flight, a feedback guidance scheme is introduced for the subsequent exoatmospheric 
stages of flight where a more simplified dynamic model permits a more analytic solution. 
Using Simplified Models 

In [2], Chandler and Smith have developed an IGM for the Saturn V vehicle. It is 
based on a flat Earth no-atmosphere model, and is further simplified with linear angle 
steering guidance. The guidance solution requires solving only a set of linear equations. 
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Ten years later, the Boeing Aerospace Company [3] adopted the linear tangent steering 
guidance as the baseline program for the Space Shuttle's PEG. Using an approximate 
gravity model, the program is extended to handle the spherical Earth case, and the solution 
is solved by an iterative algorithm. 

Perturbation Methods of Analysis 

Perturbation methods of analysis have been shown to be powerful approaches to 
spacecraft guidance design. Breakwell and Rauch [12] have used regular perturbation to 
solve a low thrust space flight problem. It is a neighboring extremal technique. A linear 
feedback control is formulated by linearizing about the reference trajectory and the solution 
is solved with a numerically determined state transition matrix. In [13], Jacobson and 
Powers have developed an explicit guidance scheme also for low thrust space flight. It is 
basically a retargeting procedure and uses an analytic solution for the inertially fixed and 
constant acceleration flight. Recently, Feeley and Speyer [14] have used regular pertur- 
bations on the expansion of the Hamilton-Jacobi-Bellman (HJB) equation, and have 
applied it to the launch vehicle guidance problem for exoatmospheric flight The approach 
requires an analytic zero order solution and quadrature evaluation. The analytic solution is 
a gain based on a flat Earth, no-atmosphere approximation, and the neglected dynamics are 
introduced as perturbations. Solution is obtained by expanding the HJB equation. In this 
method, higher order state histories are not required and higher order corrections for the 
costates are obtained by partial differentiation of the power series solution to the HJB 
equation. An alternative approach based on regular expansion of state and costates was 
also developed by Leung and Calise [15]. This approach has the advantage that on-line 
quadrature can be avoided. However, both the solution approaches of [14, 15] were later 
found to be inadequate when aerodynamic effects are included. 

1,2 Research Contributions 

The major contributions of this research are: (1) an exhaustive study and simulation 
effort which demonstrates conclusively that perturbation methods alone (both regular 
and/or singular perturbations) are inadequate for use in developing a guidance algorithm for 
the atmospheric phase of a launch vehicle trajectory, and (2) the development of a hybrid 
approach, that combines the numerical method of collocation and the analytic method of 
regular perturbation to make it suitable for real-time guidance, and superior to either method 
alone. The hybrid approach retains the desirable and complimentary features of the 
individual methods. The collocation method is further improved by providing more 
intelligent choices of the interpolation functions, which are derived from the analytically 
tractable portion of the necessary conditions for optimality. When applied to the launch 
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vehicle guidance problem, the main result is a bilinear tangent steering law for the thrust 
vector angle that can be employed for all flight phases, including the atmospheric phase of 
the trajectory. The progress reports and papers that are related to this research effort can be 
found in [15 - 25]. 

A second effort that paralleled this work under the same grant number was lead by 
Dr. Dewey Hodges, of the School of Aerospace Engineering at Georgia Tech. This work 
has been documented under a separate contractor report [26]. 

1.3 Report Organization 

Sec. 2 presents the formulation of the launch vehicle trajectory optimization 
problem, which includes the equations of motion and the vehicle aerodynamic and 
propulsion models that are based on a generic model of the ALS. The results for two 
purely analytical approaches are documented in Sec. 3. The first is a singular perturbation 
approach using an energy state approximation and a 2-state model. The second is a regular 
perturbation approach based on the zero order solution for a flat Earth no-atmosphere 
assumption. Sec. 4 details the development of a hybrid approach that employs both regular 
perturbation analysis and the method of collocation. A fourth order nonlinear system is 
treated in depth to demonstrate its application, and to compare it to solutions obtained by 
both regular perturbation analysis and purely numerical collocation methods. In Sec. 5, the 
launch vehicle guidance problem is presented using the hybrid approach. It includes the 
zero and the first order correction formulations and their solutions, and compares the 
resulting guided solution with the optimal solution obtained by the method of multiple 
shooting. Sec. 6 is the conclusions of this research and the recommendations for future 
research . 
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SYMBOLS AND ABBREVIATIONS 


Symbol 

V - airspeed 

Wy k - wind speed components in the Cijk frame 
X - heading angle 

y - flight-path angle 

m - vehicle mass 

tg - staging time (158.5s for the ALS vehicle) 

r - magnitude of radius vector measured from the Earth's center 

r e - Earth mean radius (6.378 x 10^) 

|ig - Earth gravitational constant (3.9906 x 10 14 mV 2 ) 

oo^ - Earth's rotational rate (7.27 x lO^rads* 1 ) 

h - altitude, h = r - r e 

T| - thrust throttle 

a - angle of attack, control variable in the wind frame 

P - sideslip angle, control variable in the wind frame 

M - Mach number 

c - sound speed 

c e - reference sound speed on Earth's surface (340.3ms~l) 

C D - aerodynamic drag coefficient Cp = Cp(a, M, P) 

Cl - aerodynamic lift coefficient Cl = Cl(oc, M, P) 

p - atmospheric density 

p e - reference atmospheric density on Earth's surface (1.225kgm"3) 

p - atmospheric pressure 

p e - reference atmospheric pressure on Earth's surface (101330Nnr 2 ) 

q - dynamic pressure 

Tvac - vacuum thrust 

A e - engine exit nozzle area 

S - aerodynamic reference area 

Q a - state transition matrix for the linear system A 

v - local vertical velocity component 

u - local horizontal velocity component 

0 - thrust-vector angle relative to local horizon, the control variable 

gg - gravitational acceleration on Earth's surface (gg = Pe/r e 2 ) 
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Symbols and Abbreviations (cont.) 
Symbol 


gj - small nonlinear terms (i = 1, 2) 

p x - interpolated state dynamics in the collocation formulation 

q* - interpolated costate dynamics in the collocation formulation 

Abreviations 


ALS 

HJB 

IGM 

KSC 

LEO 

PEG 

TPBVP 

MPBVP 


- Advanced Launch System 

- Hamilton- Jacobi-Bellman 

- Iterative Guidance Mode 

- Kennedy Space Center 

- Low Earth Orbit 

- Powered Explicit Guidance 

- Two-Point Boundary Value Problems 

- Multi-Point Boundary Value Problems 
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SECTION n 


PROBLEM FORMULATION 

In this section, we first formulate the optimal launch vehicle guidance problem, 
which includes the equations of motion for a point mass model of a launch vehicle that the 
subsequent analyses are applied to. The reference aerodynamic, atmospheric and propul- 
sion models are also included 

2.1 Equations of Motion 

Referring to Fig. 2.1, the point mass equations of motion for a multi-stage launch 
vehicle over a spherical, rotating Earth inside a non- stationary atmosphere are: 

V =— — cosacos P — 2 sin y + rco 2 (sin y cos 2 X- cosy sin X cos X cos x) 

m x L 

-Wj cos y sin % - Wj cosy cos % - W k sin y + 2co e [Wj(sin y cosX - 
cosysin>.cosx) + cosY(WjsinA,-Wj c cosA,)sinx] ; V(t 0 ) = V 0 

£ = {-— — costxs i n P+ X — + ^— cos 2 ytanXsinx + ro) 2 sin A. cos X. sin x 
m r 

+2co e V(cos y sin X - sin y cos X cos x) - Wj cosx + Wj sin x + 

2co e [ Wj sin X sin x + (Wj sin X - W k cos X,) cos xl} / ( V cosy) ; x(t 0 ) = Xo 

v = {— — S * n a + h. (ib. _ cos y + rco 2 (cos 2 X cosy + sin X cos X sin y cos %) 

m r z r 

+2co e V sin x cos X + Wj sin y sin x + Wj sin y cos x - W k cos y 
+2co e [ Wj (cos y cos X + sin y sin X cos x) - sin y (Wj sin X - 


W k cosX)sinx]}/ V 

/ — \ 

r-r 

O 

II 

o 

• V cosy sin x + Wj 

(() = 

rcosX 


Vcosycosx + Wj 
X — 

r 

; X(t 0 ) = X Q 

r = Vsiny-t- W k 

; r(t 0 ) = r 0 

m = f(q, r, t) ; m(t 0 ) = m 0 

; m(t ( s 1 + ) ) = 


where 
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Here y w and z w are defined in the opposite from their usual convention. 
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dW T i . 3W X ^ 


= :1II2L 

x ^ y dX 


U&Li 


3<|) T dX dh 
D (i) = qS (i) c{j } ; Y (i) = qS (i) dj } 


; x = {i, j, k} 


; L {i) = qS (i) cP 


q = pV 2 / 2 


• = nT (i) 

’ 1 l max 


; T1 € [0, 1] 


( 2 . 2 ) 


Here, an inverse-square gravitational field is assumed and }i e is the Earth's gravitational 
constant (3.9906 x 10 14 m 3 s' 2 ). A higher order harmonic model to account for the Earth's 
oblateness can be used by replacing \ijx 2 with the harmonic expression. The superscript 
(i) = { 1, 2, ... , n} indicates different stage values. The above complex model provides 
sufficient details for most trajectory analysis purposes. 

The state variables in this model are airspeed V, heading angle flight-path angle 
y, longitude 4>, latitude X, radius vector from the Earth’s center r, and vehicle mass m. The 
variables V, y, x are relative to the moving air. The wind velocity components Wj, Wj, 
are assumed to be given as functions of {(f), A., h}, where h = r - r e is the altitude and r e is 
the mean Earth radius (6.378 x lO^m). The control variables are throttle T| , angle of attack 
a and sideslip angle (3. The coefficients of drag C[>, side force Cy and lift Cl are 
functions of a, (3 and Mach number M = V/c. The fuel rate f is a function of throttle 
setting, altitude and time. The after-jettison stage mass m(t s+ ), staging time tg, are vehicle 
parameters, and are both assumed fixed here. Standard atmospheric properties such as 
density p, pressure p, and sound speed c are given functions of h. The coefficients and 
properties are given in tabular forms which are interpolated as smooth functions of the 
independent variables. 


2.2 Assumptions and Simplifications 

To simplify the analysis, the following assumptions are exercised: 

Analytic thrust expression - As mentioned in the previous section, a typical launch vehicle 
employs maximum throttle T| = 1.0 during the ascent phase. For most trajectory analysis 
purposes, thrust can be adequately modeled as 


t(0 

max 




(2-3) 


where T vac is the vacuum thrust value and A e is the engine nozzle exit area. The term A e p 
represents the back-pressure effect that causes a drop of thrust level as the engine is 
operated inside the atmosphere. 
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Constant fuel rate - For a purely rocket propulsion system the rate of fuel consumption is 
proportional to the vacuum thrust 



(2.4) 


where gg = Hg/r e 2 , and 1^ is the specific impulse, a measure of the fuel efficiency. Modem 
rocket engines have values ranging from 300s to 450s*. 

Non-rotating Earth - The Earth's rotation, to e is small (7.27 x lO'^rads' 1 ) and the term rco e 2 
which represents the transport acceleration, was neglected. The term 2o) e V which 
represents the Coriolis acceleration may reach 0.1g e at orbital speed. Here g e is the 
gravitational acceleration at the Earth’s surface. However, the vehicle reaches orbital speed 
sharply near the end of its flight phase. Therefore, the dominant effect of this term is only 
apparent for a short period of time, and setting co e = 0 does not produce any significant 
error. 

Planar motion - In actual flight, the lateral maneuver is short This magnitude is dependent 
on the launch site which is selected as close to the equator as possible so that a wide range 
of orbit inclination can be achieved. A large amount of lateral maneuver is typically not 
required and the desired flight azimuth can be achieved very early in the flight. Hence for 
simplicity, it is assumed that there is no out-of-plane motion by setting p = Wj = 0 and 
considering Cy(P = 0) = 0. These assumptions allow us to decouple the dynamics of 
airspeed, flight-path angle and altitude from those of heading angle, longitude and latitude, 
and the dynamics are reduced to those associated with motion in the vertical plane. For 
convenience, the vehicle is assumed to be launched due east on the equator, i.e. Xo = 90° 
and X 0 = 0, and <j) 0 is arbitrarily set to zero. The resultant system is a 4-state model: 


V = 


T^cosa-D^ 

m(t) 


Lip 

-^-f-siny- W, cosy - W k siny 
r 


; V(t Q ) = V 0 


y = { 


T (i) sina + L (i) 
m(t) 


: Vcosy + W; 

<i> = 

r 


r = Vsiny + W k 



V 2 • • 1 

— )cosy + Wi siny- W k cosy}— ; y(t Q ) = y 0 

; <t>(t 0 ) = 0 

» r(t 0 ) = h 0 + *e (2.5) 


where 


* For comparison, specific impulse of tubojet engine is over 5000s. 
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n-1 


m(t) = 


fm 0 -k (1) (t-t 0 ) ;t 0 <t<4 1} 

i< j) - k<J)(t - #) ; < t < t< j+1) ; t< n) = t f ; j = 1, 


k (i) =T® /(gel®) 


( 2 . 6 ) 


The initial conditions chosen for this problem represent the vehicle states following 
a vertical launch and clearing of the launch tower. The terminal constraints represent direct 
injection at the perigee of an 80nm x 150nm elliptical transfer orbit. 

V 0 = 64.49m / s ; y 0 = 89.5° ; h 0 = 400m ; t 0 = 15s 

V f = 7858.2m /s ; y { = 0° ;h f = 148160m (2.7) 

The objective is to minimize the final time, which is equivalent to minimizing the fuel 
consumption for this formulation. Since there is no constraint on <|>f T the <{> dynamics in Eq. 
2.5 are ignorable and can be deleted from the analysis. Also, the optimization must be 
performed subject to the constraints q < q m ax and I otq | < (otq) max . 


2.3 Aerodynamic Model and Launch Vehicle Configuration 

The aerodynamic model (cf. Figs. 2.3 - 2.8) is obtained from [27]. It corresponds 
to a generic model of a heavy-lift capacity 2-stage launch vehicle based on a CFD analysis. 
The vehicle has an asymmetric configuration as shown in Fig. 2.2 with the booster 
mounted atop the main body. The booster produces a shadowing effect above supersonic 
speeds during the first stage flight. This shadow effect reduces the Cj) at positive angle of 
attack and the Cq exhibits a nonconvex behavior (cf. Fig. 2.6) in a above Mach 1.3. 
Other than this behavior, Cu(a) and Cl(oc) are nearly parabolic and linear respectively at all 
Mach numbers. 
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Due to the nonconvexity in C]> the Hamiltonian also becomes nonconvex. The 
control is expected to jump as a switches from a lower value to higher value when the two 
peaks (using Maximum Principle) of the Hamiltonian become equal as time progresses. 
The phenomenon is displayed in Fig. 2.9. The study documented in [22] has shown that 
the hodograph can be convexized when bank angle is used as a second control variable, 
and the angle simply switches from 0 to 7t to make use of the lower Cd at small positive a. 
It has also shown that the effect of a chattering control of the first kind [27], if exists, will 
be small and that a chattering arc is not expected. This hypothesis is consolidated by the 
numerical analysis here, where no high frequency control activity is observed within the 
nonconvex region. 




in Control due to Nonconvex Hamiltonian. 


Table 2.1. ALS Vehicle Physical Data 



1st- stage 

2nd-stage 

mo(to); m s (ts + ) 

1,523,400kg (15s) 

546,600kg (158.5s) 

T 

1 vac 

25,813,000N 

7,744,000N 

^sp 

430s 

430s 

S 

131.34m 2 

65.67m 2 

A 

37.51m 2 

11.25m 2 

Qmax 

40698 .2Nm -2 

nil 

(Otq)max 

167,580degNm -2 

nil 


Since sideslip is not considered, the aerodynamic coefficients can be interpolated as bicubic 
splines [28] in a and M. The interpolation scheme provides up to second order continuous 
derivatives. Other physical parameters of the ALS vehicle arc given in Table 2.1. 

2.4 Atmospheric Model 

The atmospheric model is based on the 1975 U. S. Standard Atmosphere [29]. 
Profiles of normalized density, pressure and sound speed with respect to their reference 
values at the Earth's surface (p e = 1.225kgm"3, p e = 101330Nm -2 , c e = 340.3ms -1 ) are 
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given in Fig. 2.10. To investigate the effect of wind shear, a mean winter wind profile 
over Kennedy Space Center (KSC) is used to model the non-stationary atmosphere. The 
profile is shown in Fig. 2.1 1. It indicates a head-on wind for vehicle launched due east, 
and the vertical and horizontal (north) wind speed components are assumed to be zero. 




Altitude (m) 

Figure 2.11. KSC Mean Wind Profile. 
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SECTION III 


ANALYTICAL APPROACHES 

Two analytical approaches are presented with the objective of simplifying the 
optimal guidance problem described in Sec. 2. The analytical and numerical results are 
summarized in this chapter. The analysis results are for: (1) a singular perturbation 
formulation, and ( 2 ) a regular perturbation formulation. 

3.1 Singular Perturbations 

Singular Perturbation theory is related to the study of a reduced solution of 
singularly perturbed systems of O. D. E's and the construction of a matched asymptotic 
series representation of the exact solution. For example, consider the following initial 
value problem 

dx 

— = f(x,y,t) ;x(e, 0 ) = x o 

dt 

£- 7 - = g(x, y, t) ; y(e, 0 ) = y 0 (3. 1) 

dt 

where x and y are scalar functions and £ > 0 is a scalar parameter. Setting £ to zero, we 
have the reduced system. Generally the reduced solution will not satisfy initial conditions 
on y, and the initial behavior of the reduced solution will be quite different from that of the 
exact one. This loss of boundary conditions on y (meaning that the reduced solution does 
not provide a uniformly valid approximation for y) is a characteristic of singular 
perturbation problem formulations. Basically, the system is separated into the slow 
variables of x and the fast variables of y. The reduction of higher order problems into 
lower order ones and the separation of numerically stiff parts by using different time scales 
are the main advantages of the method. Applications of the method are detailed in [30, 31]. 
a) Energy state approximation 

The energy state approximation is the most widely used approximation in aircraft 
performance optimization, and sometimes referred to as energy management. It has been 
applied to minimum time-to-climb, minimum fuel-to-climb and minimum time intercept 
problems. First we replace the velocity with the mass specific energy 


* A third analytical attempt using matched asymptotic methods is documented in [25]. 
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(3.2) 


E = V 2 /2-|i e /r 


as the state variable. Differentiating Eq. 3.2 and using Eqs. 2.5, 2.6 leads to the system 

• T^cosa-D^_, 

E = V 

m(t) 


r = Vsiny 



sin a + 
m(t) 


~(h~— )cosy}/V 
r r 


(3.3) 


where V - ^2(E + q e / r ). At the moment, the wind shear effects are not considered. In 
earlier studies on supersonic aircraft [30, 31], specific energy and mass are regarded as 
slow variables and altitude and flight-path angle are treated as fast variables. So to put Eq. 
3.3 into the singular perturbation form, we artificially introduce a bookkeeping parameter 8 
into Eq. 3.4 as follows: 

^ T^cosa-D^ __ 

E = V 

m(t) 


er = Vsiny 


ttO) 
ey = { 


sina + li 1 ^ 
m(t) 


-(%-— )cosy}/V 
r r 


(3.4) 


The performance objective is to minimize tf. 

The necessary conditions are formulated by first moving e to the right hand side of 
the differential equations, and define the Hamiltonian as 


^ T^cosa-D^\ Ifi, • f T^sina + l/^ 

H = A E V +— Vsiny + —{ 

m(t) e e m(t) 

M V 2 

-(tt )cosy} / V + constraints 

r r 


(3.5) 


The costate dynamics satisfy: 


Ae 


9H 

3E 




9H 

dy 


(3-6) 
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(3.7) 


Now introduce the transformations - ^E» e ^r ^r» e ^7 ^"Y which results in 




3H 

BE 




m 

dy 


Note that X E is a slow variable and that ^ ^ are fast variables. The optimality condition 
is given by 

9H/3a = 0 ( 3 - 8 > 


In the reduced problem (e = 0) r and y are treated as control-like variables, which is a 
consequence of setting e = 0 in Eq. 3.7. The transformed costates A, r , Xy (when 
substituted in Eq. 3.6) can be interpreted as Lagrange’s multipliers used to enforce the 
constraints that result from setting e = 0 in Eq. 3.5. 

Reduced (outer) solution 

The reduced or outer solution corresponds to the solution of Eqs. 3.4, 3.7 and 3.8 
when e is set to zero. The condition dH/9r = 0 (which results from setting £ = 0 in Eq. 
3.7) is a first order necessary condition for a minimum of the Hamiltonian in Eq. 3.5 (we 
are minimizing the final time). Since the costate Xg may be interpreted as 0tf/3E(t o ), it 
follows that in the reduced problem, X E < 0. Hence a stronger statement for this optimality 
condition may be written as 


♦ 

r 


max 

r 


T (,) cosa - D (l) 
m(t) 



E,m 


(3.9) 


subject to the conditions: 


y = 0 

T (i) sina + L (i) 


0 = 


m(t) 


|! e V\ 

-(rr )cosy 

r z r 


q < 40698. 2Nm -2 
M < 2924.82radNm' 2 


(3.10) 


The last two conditions are the dynamic pressure and aerodynamic load constraints. 
Starting at an initial energy level and initial mass, a one-dimensional search in altitude is 
performed. The energy level is then increased and the corresponding change in mass is 
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V (m/s) 

Figure 3.2. Angle of Attack Profile along the Reduced Solution. 







Since the optimal solution also exhibits a large value of flight-path angle (incon- 
sistent with the reduced solution approximation), another calculation scheme is used to 
estimated as 

Am = -k (i) (AE / E*) (3.11) 


where the superscript '*' denotes evaluation on the reduced solution. Hence by sweeping 
through all the energy levels of interest, a reduced feedback guidance law that defines the 
best altitude profile is obtained. Figs. 3.1 and 3.2 show the results for the reduced 
problem when the optimization in Eq. 3.9 is carried out for the first-stage flight. The initial 
conditions in E and m are chosen along a reference optimal trajectory. The solutions at low 
energy levels result in very large values of angle of attack (> 20°) that are well beyond the 
given aerodynamic model range and therefore should not be considered feasible. The 
reduced solution is unrealistic in that the vehicle stays on the <xq constraint up to an energy 

level of -6.09 x 10 7 Jkg' 1 . 

Since the optimal solution also exhibits a large value of flight-path angle (incon- 
sistent with the reduced solution approximation), another calculation scheme is used to 
estimate a non-zero flight-path angle and to include the effect of a non-zero flight-path angle 
in the reduced solution. Assuming the vehicle is already on the reduced solution and is to 
follow the trajectory, the change in altitude along the reduced solution gives an estimate of 
the flight-path angle according to 


sin Y e = 


Ah 

(AE / E)V 


* 


( 3 . 12 ) 


By perturbing the energy level from E to E + AE, we have 

Ah* = h* (E + AE) - h* (E) (3. 1 3) 

and a central difference scheme is used to estimate y. Then the solution of Eqs. 3.9 and 
3. 10 is recalculated with y = 0 replaced with y = y t . The results are given in Figs. 3.3 and 
3.4. The inclusion of y t gives a slightly lower value of a, and both angle profiles behave 
reasonably. However, this calculation scheme becomes numerically unstable once the 
vehicle left the aq constraint boundary. 
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b) Two-state model 

Since the energy state approximation does not produce a solution that resembles a 
reasonable flight trajectory, a more accurate model is employed. The new reduced-order 
model corresponds to a 2- state approximation: 


E = 


'p(i) _ p)(0 

m(t) 


V 


r = Vsiny 


T (l) a + K{ l} a u e V 2 

ey = { (% )cosy}/ V 

m(t) r r 


(3.14) 


where only the flight-path angle is assumed fast. To make Eq. 2.5 analytically tractable, 
we adopt the assumptions that the induced drag due to a is negligible and lift is linearly 
proportional to a (L® = K^a). The necessary conditions for optimality of flight-path 
angle and angle of attack on the reduced solution are: 


BH 

dy 


= >, r Vsiny 


31/2 

;X T <0 

singular 

;X T =0 

-k/2 

\X T >0 


a = 


m(t) 
T (i) + 


(%-— )cosy 
r z r 


(3.15) 


In [23] it is shown that the velocity hodograph for the 3-state reduced model (including 
mass) is nonconvex, and that at A. r = 0 the optimal solution chatters between y = ±k/2. The 
interpretation here is that when the altitude reaches its optimum value (for the current 
energy and mass), then a chattering solution is able to maintain the optimum altitude rate 
while maximizing the ratio of the mass rate to energy rate. Therefore this formulation is 
totally inappropriate for the analysis of energy climb in that it produces a reduced solution 
made up of vertical climbs and dives, connected by chattering arcs, 
c) Manifold solution and eigenvalue analysis 

The fundamental problem inherent in treating launch vehicle dynamics by energy 
state approximation relates to the constraints on the y and h dynamics. They are fast in 
comparison to energy and mass dynamics and without taking into account the dependency 
on the singular perturbation parameter e. For instance, the constraint on altitude dynamics 
implies y = 0 along the reduced solution, which is an extremely crude approximation for the 
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Time (s) 

Figure 3.5. Flight-path Angle Profile for Various Reference Trajectories. 



Time (s) 

Figure 3.6. Angle of Attack Profile for Various Reference Trajectories. 
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launch vehicle case. This problem can be alleviated by using a slow manifold solution [32] 
in place of the reduced solution, which amounts to solving the exact problem with the initial 
flight-path angle chosen to suppress any fast motion that may be present in the solution. A 
separate boundary layer analysis could then perform to take into account the actual initial 
condition on y (cf. Fig. 3.7). This approach has also been carried out, however it is found 
that the assumptions regarding the separation of dynamics worsen above supersonic speed, 
and the reduced-order model approximation deteriorates. This hypothesis is consolidated 
by the eigenvalues investigation described below. 

Computation of the equilibrium manifold corresponds to determining the initial 
condition on y so that rapid transients in y and Xy are absent in the exact solution. First a 
sweep of the initial condition in y about a nominal value of y Q = 89.5° is performed, and the 
exact dynamics of the states and costates (with the control eliminated using the optimality 
condition) are numerically integrated. This allows us to identify the equilibrium manifold 
by visual inspection for the absence of fast transients in y and Xy. The closer the actual 
initial condition for the fast variable lies to the manifold, the more accurate the subsequent 
boundary layer correction in y becomes. Figs 3.5 and 3.6 demonstrate that the manifold is 
estimated to be at y 0 = 75°, where it can be seen that there is no apparent boundary-layer- 
like behavior in the fast variable y and the control a. 



Figure 3.7. Typical Boundary Layer Characteristics. 

To shed insight on the separation phenomenon of the fast and slow dynamics of the 
launch vehicle problem, an eigenvalue test is carried out By linearizing the dynamics of E, 
r, y, A.£, Xf, Xy about the equilibrium manifold, the eigenvalues of the linearized system are 
obtained, and the relative magnitudes of the real part of the eigenvalues provide information 
about the separation possibility of the dynamics. A Hamiltonian matrix appears in the 
linearized system whose eigenvalues characterize the full order system of dynamics (states 
and costates) in the vicinity of the equilibrium manifold. 

Eigenvalues calculated at discrete points along the trajectory are shown in Figs. 3.8 
and 3.9 (only those in the right half s-plane are shown). At the beginning part of the trajec- 
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tory (t < 50s), the results clearly show a separation configuration of 2 slow and 1 fast state 
(and costate) variables. All the eigenvalues are real. The relative magnitude is separated by 
a factor of up to 4 in this interval (cf. Fig. 3.9). As the energy level increases, two of the 



Time (s) 


Figure 3.8. Eigenvalue Analysis along the Reference Trajectory of y Q = 75°. 



Time (s) 

Figure 3.9. Eigenvalue Separation by Relative Magnitude. 
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eigenvalues join to form a complex conjugate pair, the real part of which is an order of 
magnitude larger than the third (real) root This suggests a decomposition of 1 slow and 2 
fast state variables. An eigenvector analysis indicates that the fast state variable at low 
energy levels indeed corresponds to the flight-path angle, whereas at high energy levels 
specific energy is the only slow state variable. Altitude, which was a slow variable at low 
energy levels, rapidly transitions to being a fast variable at approximately t = 50s as shown 
in Figs 3.8 and 3.9. 

A nonlinear feedback control solution for angle of attack, based on a boundary layer 
correction for the flight-path angle dynamics, can be formulated as follows: 

H = — X, mo k^^ + A,£ 0 E(E 0 ,h 0 ,m 0 ,cx) + X. r0 V 0 sinY + X,yY(E o ,h o ,m 0 ,(x) — 0 

H a = 0 (3-16) 

where itiq, E q , h Q , X, Eo , Xj- 0 are treated as slow variables* in the manifold solution, 
and are constant in the boundary layer analysis. The manifold solution is stored as a 
function of energy, and the boundary layer problem defined in Eq. 3.16 is solved at each 
control update to form a guided solution. Note that there are two equations for the two 
unknowns in a and Xy. The guided solution using the pre-computed slow manifold 
(chosen for y 0 = 75° in Fig. 3.6) with an on-line boundary layer correction is plotted in 
Fig. 3.10. The optimal solution approaches the manifold solution. However the guided 
solution is first attracted to the manifold, and then diverges at about t = 25s. This correlates 
almost exactly with the transition that takes place in the eigenvalue associated with the 
altitude state in Fig. 3.8. That is, the role of altitude variable has changed, but the 
boundary layer analysis has treated the altitude variable as slow (constant to zero-order in 
e). This explains the failure of the manifold approach for this problem. 

Recalling the previous energy state approximation formulation, even though eigen- 
values analysis clearly indicates the existence of a two-time scale behavior, the poor 
performance of the zero order reduced solution is attributed to the large value of 
longitudinal load factor inherent in the launch vehicle problem. The value of this non- 
dimensional variable along the reduced solution is plotted in Fig. 3.10. In comparison with 
a subsonic transport aircraft with a load factor of 0.1, the launch vehicle averages above 3 
in this case and therefore a zero or even a first order solution is not expected to provide any 
reasonable approximation. 


* Here we treat m as a state variable, and it is not eliminated by Eq. 2.6. 
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3.2 Regular Perturbations 

The unsuccessful attempt by singular perturbation analysis led to consideration of 
another analytic approach that has been used repeatedly on low thrust spaceflight problems, 
the regular perturbation analysis. In this section, the general regular perturbation 
formulation for optimal control problems is discussed. An extension over earlier 
formulations is that higher order corrections for the free final time are made explicidy in the 
formulation developed here. Then an analytic zero order solution based on the maximum 
horizontal speed transfer problem in a constant gravity field and in vacuum [33] is extended 
to a mass-varying multi-stage rocket. This is then followed by an attempt to compute a first 
order correction to account for a central gravitational field, spherical Earth and all the 
atmospheric effects. 

al Re gular perturbations in optimal control 

The optimal control problem formulation consider here is to maximize a perfor- 
mance index which is a function of the terminal states and time, subject to dynamic 
constraints: 


J = 


max 

u 


{<t>(x,t)} 


(3.17) 


x = f(x,u,t) + eg(x,u,t) ; x(t Q ) = x 0 ; t e [t 0 ,tf] 


(3.18) 


and the terminal time constraints \jq(x(tf)) = 0, i = 1, ... , p ^ n. In Eq. 3.18, x is an n- 
dimensional state vector and u is an m-dimensional control vector. In applications, the 
expansion parameter e is sometimes artificially inserted to signify the presence of small 
nonlinear effects, and used as a bookkeeping parameter for the regular expansion analysis. 
The Hamiltonian and transversality condition are given by: 


H = X T {f + eg} 


; H(t f ) = -d> t 


O = <() + v 1 y 


(3.19) 


The costate equations and associated boundary conditions are: 


X = 


" H x 


; X(tf) = 4>x 


tf 


(3.20) 


where the subscript is used to denote partial differentiation. In the absence of control 
constraints, the optimal control satisfies 

H u =X T {f u +eg u } = 0 (3.21) 
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assuming that H uu > 0. 

In the above final time is free. Thus, we introduce a new independent variable x = 
(t - to)/T where T = tf - ^ and rewrite the necessary conditions of Eqs. 3.18 - 3.20 in the 
following equivalent form: 

*' = H X T ; x(t = 0) = x 0 ; \jt(x(x = 1» = 0 (3. 22) 


A/ = -H X T 


; Mx = 1) = <D X 


x = 1 


(3.23) 


r = 0 


(3.24) 


H = X. T {f(x,u,xT + t 0 ) + £g(x,u,xT + t 0 )} ; H(x = 1) = -O t 


x = 1 


(3.25) 


where (-)’ denotes d(-)/dx. In a regular perturbation analysis, the objective is to approxi- 
mate the solution to Eqs. 3.22 - 3.25 by an asymptotic series in x, A., u and T as follows: 

x = xq + exj + e 2 X 2 + ... 

A, = Xq + eA,j + £ 2 A,2 + ... 

U = Uq + £Uj + £ 2 U2 + ... 

T = To + eTj + £ 2 T 2 + ... (3.26) 


Assume the functions f, g, <)>, \j/ have piecewise continuous derivatives up to order at least 
K+l where K is the order of approximation. Using the Taylor series formula, a finite 
series approximation is constructed according to 


F(o 0 +ia k £ k ) = F(a 0 )+E^ 
k=i da 


2 ( dF 
<*i+e {— 
<*0 da 


1 d 2 F 


a 2 + 9 

o 0 1 2! da 2 


Po 


a 2 }+ 


(3.27) 


where a = {x, X,, v, u, T}. Substituting the series representation for each of the variables 
in Eqs. 3.22 - 3.25 and equating like powers in £, we obtain the zero order and higher 
order necessary conditions. To zero order we have: 

9x 0 / at = aH 0 / ax. 0 ; x 0 (t o ) = x 0 ; \|/(x 0 (t o + T 0 )) = 0 


ax. 0 / at = -aH 0 /ax 0 ; x. 0 (t = t 0 +t 0 ) = <D(x 0 ,i)/ax 0 


1 = l o + T 0 
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8H 0 / du 0 = 0 


H 0 = ^ T f(x 0 ,u 0 ,t) 


; H 0 (t = t 0 +T 0 ) = -O(x 0 ,t)/at 


t = t G + T 0 


(3.28) 


In Eq. 3.28, the new independent variable t = tTq has been introduced, where it should be 
noted that in the zero order problem T = Tq. 

For the higher order problems, they are governed by a set of nonhomogeneous 
linear O. D. E's. with the form of 


d 

’ x k' 


dt 

X k _ 



A|i(xo.Xo.T 0 ) A 12( x 0 

A 2l( x 0’^0’ T 0) A 22( x 0 


A 0 .T 0 ) x k Tk_ C l( x 0’^0’ T 0) 
>^0’To). _^k. Tq _C2( x 0’^0’T())_ 


Plk( x oA 0 ,T 0 , — ,Xk-i,A.k-i,Tk_i) 
P 2k( x 0’^0’ T 0’ — > x k-l’kk-l’ T k-l) 




(3.29) 


where 

An = f x -fuKfJX^rkfJWx 
Ai2=-f u t(f u T >-)ur 1 f„ T 
a 2 i = -(fjx) x 

a 22 = -H +(fx T x) u [(fXr‘f u T 

C, = f + (t - t 0 ){f; - f 0 ((f g T X) u ]-‘(f„ T ).)- t } 

C 2 =-f^-(;-t 0 ){(f x T X);-(fx T X) u [(fT>.) li r 1 (fu T >.) ; } (3.30) 

and for k = 1 : 

Pl^g-WuX]-^ 

P 21 = -gIX + (f x T X) u [(fJX) u r 1 gJX (3.31) 

All the matrices in Eq. 3.30 are evaluated at the zero order solution values. To complete the 
necessary conditions, it is also required to expand the boundary conditions and the 
transversality condition in Eq. 3.28. Note that Eq. 3.29 explicitly shows the effect of 
higher order corrections to the final time, T. If the solution process is terminated at say. 
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k = 1, then a real-time sampled data implementation of the control solution would be 
constructed as follows. For the original system in Eq. 3.18, an expression for the optimal 
control is obtained as function of x and X from the optimality condition in Eq. 3.21. Then, 

treating the present state as the initial state, a first order approximation is obtained by using 
^o(to) + ^l(to) as an approximation for A,(t 0 ) to compute the control, where A.g(t 0 ) and 
^-l(to) are obtained from the solutions of the zero and the first order necessary conditions. 
This process is repeated at the next control update time by regarding the value of the state as 
the new initial state. Therefore, it is necessary to repeat the zero and first order solutions in 
updating the estimate of the costate variable. 

The non-homogeneous linear ordinary differential equations in Eqs. 3.29 - 3.31 
may be expressed in terms of a convolution by first obtaining a state transition matrix. The 

A 

state transition matrix i2^(t,t 0 ) is merely the partial derivative of the zero order solution at 
t with respect to the initial conditions xg^) and Xg(to), hence it is easily computed given 

an analytic zero order solution. In Appendix A it is shown that the result can be expressed 
in the following form 


x k(0 

Mt) 




X k (t 0 ) 

Ak(to). 


+ ^o Q A(t.t)| T ^ 


Qd)' 

c 2 (x) 





x k(t 0 ) 

AkOo) 


+ T k 


t~tr 


T 0 


x o(0 


„ rp lk (T)- 

+J t l Q A (t,x) 1 d 
° L P 2k( x )J 


it (3.32) 


Using the above expression at t = Tg along with the expansions of the boundary condi- 
tions, we can solve for X k (t 0 ), v k and T k from a set of linear algebraic equations. Thus the 
major part of the computation for the first order term lies in the quadrature that must be 
performed in Eq. 3.32. In a discrete time implementation, if the current state is regarded as 
the initial state then xj c (t Q ) = 0 in Eq. 3.32 since xg(to) satisfies the initial condition on the 
state variable. Since zero order solution changes at each update of the initial state, it is 
necessary to repeat the quadrature at each update for the higher order corrections. 
Alternatively, we can fix the zero order solution and treat x k (to) as the deviation between 
the current state and the zero order solution computed for the original epoch time, but 
evaluated at the present time. In this form it would be possible to pre-compute the 
quadrature and store it as a function of a monotonic variable along the trajectory. Thus the 
real-time process of solving the zero order problem and the quadrature can be avoided. 

The case of discontinuous dynamics, such as might arise in a multi-stage launch 
vehicle, can be handled by a simple modification of Eq. 3.32. For example, in a two-stage 
representation we would have 
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(3.33) 


x k (t) 

M t) 


= n ( 2>(it s ){ Q ( ^(t s ,t 0 ) 
jJsQ^Ct,!) 
j‘n ( A 2) a,x) 


pf?(t) 


lk 
>(2 

•ffw. 


x k^o) 

[_A,k(to)J 

dx } + T k 


+ T k 


t-t c 


tc-tf 




.^(ts)] 


To 

|dx ; t > t s 


4 2 W 

|> ( 0 2 ) <' t>J 


The superscripts (1), (2) denote the expressions for different sets of dynamics and tg is the 

interior point where discontinuity occurs, 
b) Launch vehicle application 

The performance objective is to maximize -tf (ie. minimize final time) subject to the 
terminal conditions V(t f ) = 7858.2ms' 1 , y(t f ) = 0, h(t f ) = 148160m, open <}>(t f ). These 
conditions correspond to direct injection at the perigee of an 80nm x 150nm elliptical 
transfer orbit. First, it is necessary to derive a closed form, zero order solution which 
should be simple, but accurate enough such that the neglected dynamics can be corrected in 
a first order term. 

Assuming that the dominant forces on the launch vehicle are thrust and gravity, an 
attempt is made to treat the atmospheric effects as a perturbation effect To further simplify 
the problem, spherical Earth effects are also considered as perturbations (these effects are 
only apparent when the vehicle reaches orbital speed near the end of the flight). The result 
is similar to the maximum horizontal speed transfer problem in [33] for a flat Earth no- 
atmosphere situation. The differences here are that the mass of the vehicle is varying, the 
dynamics are discontinuous and the terminal boundary conditions are specified at an 
unknown final time. We now recast the dynamics of Eq. 2.5 in a regular perturbation 
format as in Eq. 3.18, in accordance with the above desired approximations: 


T^ sin 0 

v= 7? — -g e +g 


m 


®-k«t 


f -pA^ sin 0 - D (l) sin y + L W cosy 
m«-k^ 


V 


. „ , U 
ge "7 r 7 


2\ 


> V(t 0 ) — v o ’ * T 2 


u = 


T«cose 


m 


(i)_ k (i) t 


+ e 


— pA^cos0 — D^ cosy- sin y uv 


(i), 


m 


(i)_ k (0 t 


; u(t 0 ) = u 0 


r = v 


> r(t 0 ) — r o 


(3.34) 
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where 


= m 0 + k (1) t 0 ; m^ 2) = m s + k (2) t s 


v = V sin y 


; u = Vcosy 


; 0 = a + y 


(3.35) 


Here e has been artificially introduced as an arbitrary bookkeeping parameter. The 
dynamics are expressed in a rectangular coordinate system to facilitate the closed form 
derivation of the zero order solution. The state variables v and u are the local vertical and 
horizontal velocity components. The control variable is 0, the thmst- vector angle measured 
from the local horizon. 

The necessary conditions of optimality for the above formulation are: 

i T — x r +e^-x v M_x 11 Mj 


i 


u 


=e Kfr 







0g2^ 

dr 


0 = (^, v cos0 - A. u sin 0) 


tO) 

vac 






+ A,,, 



0 = {^. V v + + ^i-r} 


-1 


(3.36) 


where the last two are the optimality and the transversality conditions respectively, and 
_ -pA^sin0-D^ sin Y + li 1 ) cosy , _ fi e u 2 

Si (!) TTfr + 8e T + 

nr ' - k v ; t r r 


62 = 


-pA^ cos0-D^ cosy ~L^ sin y uv 

m (i) -k (i) t 7 


(3.37) 


Zero order solution 

Setting e = 0, the costate solutions and the optimal control are given as follows 
(with some license taken with the zero order time notation): 

^vo(f) = c v0 — c r0* » ^uO(f) = c u0 > ^rO(f) = c r0 
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tan(0 o (t)) = p = qt ; p = c v0 / c u0 


5 q — c rO I c uO 


(3.38) 


The control satisfies a linear tangent law. Substituting Eq. 3.38 into Eq. 3.36, the state 
equations can be integrated in closed form. The solution that relates the states at t > t§ to the 
initial conditions is presented below and for t < tj, the terms involving variables with 
superscript (1) would simply be deleted. 


( 1 ) 


vo(t) = v 0 -g e (t-t 0 ) + T ^G(m (1) ,k (1) ,x) 


T = t s 
X = t n 


r(2) 


^-G(m (2) ,k (2) ,x) 


r( 2) 


X = t 
X = U 


r(l) 


u 0 (t) = u 0 + T ^-F(m (1) ,k (1) ,x> 


T = t s 
x = t r 


y(2) 


7§£-F(m (2) ,k (2 \x) 


r(2) 


X = t 
X = U 


( 1 ) 


1 | v x / 

r 0 (t) = r 0 + v 0 (t - t Q ) - -g e (t - 1 0 ) 2 - (t - t 0 )-^G(m (1) ,k (1) ,t 0 ) - 


f T (2) T (l) 

(t - t s > -4g-G(m (2) ,k (2) ,t s ) - -^-G(m (1) ,k (1) ,t s > 

>. - 


T d) 

^ffrK(m (1) ,k (1) ,x)| 
qk' ; 


X ts +i^K(m (2) ,k (2) ,x)| 


x = t 0 qk 


f + 

x = t 

X = to 


; t > t s (3.39) 


where 


F(m (i) ,k (i) ,x) = 


-sinh 1 [tan(0 o (x)-q)] 4 1 . A qm (l) - pk (l) 

7 ,tanq- A ,A- (i) 




+ A Z 


G(m (i) ,k (i) ,x) = -AF(m (i) ,k (i) ,x) - sinh -1 [tan(0o(x))] 

K(m (i) ,k (i) ,x) = -sec(0 o (x)) - [tan(0 o (x)) + A]G(m (i) ,k (i) ,x) (3.40) 


To solve for the solution, Eqs. 3.38 - 3.40 are evaluated at the zero order final time tfg = to 
+ tg + Tq where Tq represents the zero order, second stage, open flight time, and used to 
enforce the zero order expansion of the terminal boundary conditions and the transversality 
condition given below: 


v (tfo) = v f ; u(tfo) = Uf ; r(tfo) - h f + r e 


, , T< 2 >sin0 x . T^cosO 

^ v ° ( m (2) _ k (2) t ge) + m (2) _ k (2) t + X r0 v 0 


|tfO 


= 1 


(3.41) 
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There is a total of four unknowns c^, c u0 , c^, T 0 to be evaluated by the four conditions in 
Eq. 3.41. 

First order solution 

Using Eqs. 3.29 to 3.31, the first order correction dynamics for the launch vehicle 
problem become 


d_ 

dt 
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r l 


A, v i 


^ul 


_A. r i_ 



0 0 0 ag 


,(i) 

*15 


0 0 0 
1 0 0 
0 0 0 
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+ 

0 
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0 
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pg’C) 


(3.42) 


with Vl (to) = u^) = r 1 (t 0 ) = vjCtfo) = u^) = r^tffl) = 0, where 


•pO) 

vac / 

v 


cos 2 0 


a 14 = — 

nr 1 ' _ j c (0 t x, v0 sin 0 q + A. u o cos 0 q 
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(3.43) 
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All the variables are evaluated along the zero order solution. Since the first stage flight time 
is assumed to be fixed, T = tf - tj. Consequently, Tj = 0 for the dynamics describing t < t s , 
and the second term in Eq. 3.42 is discarded for the correction dynamics corresponding to 
this time interval. In this example, the state transition matrix has a structure of 


£2^(12. tl) = 


1 

0 

0 


(VjO) 

“l5 

1 

O'O 

3 

0 

1 

0 

“24 

m (i) 

“25 

“26 

l 2 -t l 

0 

1 

m (l) 

c0 34 

m (i) 

“35 

“36 

0 

0 

0 

1 

0 

l l -t 2 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 


(3.44) 


Complete expression of the co's are given in Appendix B, and the expansion of the trans- 
versality condition for the first order case is 


0 = 


! - 8= + ao>- + + g2o) + 

^(2)^(2) 

T 1 — wr m 2 ^v0 sine 0 + ^ U 0 COS0 O ) 

(m w - k (Z) t r 


t = tf0 


(3.45) 


From Eq. 3.33, the first order variables at tfo are related to their initial values at ^ by 


x l(tfo) 

Ai( t fo)J 


= ^A^(tfO’t S ) 


*P(lfO> 




Xl(t 0 ) 

LXiOo). 


“1 t 0 




+ t J 0 £2 ( A , (tf0.X)| 

to 


p(2 ) 

Ml 

P (2) , 
“21 J 


o 

Idt 


p(l) 
Ml 
,p(D 
L/21 . 


Idt 




(3.46) 


where xj = [vp Up rj }, = {^vi» ^ u l> ^rl}> ^11 = tPl* P2> ^21 = tP4> P5’ P6^- 

Substituting Eq. 3.46 into Eq. 3.45 and using the boundary conditions defined in Eq. 
3.42, the unknown costate and final time corrections ^vi(to), A. ul (t 0 ), )tri(t 0 ), Tj can be 
found from a set of linear algebraic equations. 

Figures 3. 12 to 3.18 give the zero and first order results for a no-aerodynamic force 
case (obtained by setting the reference area S = 0). The optimal solution obtained from a 
multiple shooting code [10] is also included for comparison. As far as spherical Earth and 
back-pressure effects are concerned, the regular perturbation approach produces very 
accurate results, especially in the state histories. Next the aerodynamic effect is included. 
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the resulting angle of attack profile is shown in Fig. 3.19. No reasonable first order 
solution is found at low altitudes in the region of high dynamic pressure and aerodynamic 
forces. The first order solution over-corrects the zero order result and gives a very large 
value of angle of attack that is not considered feasible. The conclusion that is drawn from 
these results has been that the aerodynamic forces are simply too large to be ignored in the 
zero order solution. Figure 3.20 show the ratios of the aerodynamic forces to the thrust 
components along the optimal solution. The magnitude of lift to thrust ratio reaches almost 
40% over some time interval during the first stage flight and indicates that a significant 
amount of aerodynamics effects exist in the ALS vehicle. 



Figure 3.12. Regular Perturbation Results in v with 

Spherical Earth and Back-pressure Effects. 
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Figure 3. 13. Regular Perturbation Results in u with 

Spherical Earth and Back-pressure Effects. 



Figure 3.14. Regular Perturbation Results in h with 

Spherical Earth and Back-pressure Effects. 
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Figure 3. 15. Regular Perturbation Results in \ with 

Spherical Earth and Back-pressure Effects. 
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Figure 3. 16. Regular Perturbation Results in with 

Spherical Earth and Back-pressure Effects. 
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Figure 3. 17. Regular Perturbation Results in ^ with 

Spherical Earth and Back-pressure Effects. 
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Figure 3.18. Regular Perturbation Results in a with 

Spherical Earth and Back-pressure Effects. 
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SECTION IV 


A HYBRID COLLOCATION/REGULAR 
PERTURBATION ANALYSIS 

This chapter develops a solution approach for nonlinear optimization problems that 
seeks to combine the desirable features of analytical methods which are based on the use of 
simplified models, and numerical methods that use elementary interpolation functions and 
finite elements to represent the solution. The approach is developed for a combination of 
regular perturbation analysis and collocation technique. A simple fourth order nonlinear 
system is used to illustrate the conceptual approach for several possible levels of 
approximation. 

4.1 Introduction 

Among the proposed analytical approaches for real time guidance in Chapter 3, the 
analysis by regular perturbation expansion of the solution is most appealing. However, 
crucial to the success of the method is that the optimal solution is reasonably approximated 
by the zero order solution, so that the addition of first or higher order corrections to the 
series solution (which usually is not convergent) results in an improvement in accuracy. 
The approach has had great success when applied to systems with small nonlinear terms 
[34, 35] so that the zero order problem is linear. Also, in certain applications a state 
transition matrix may be determined for the first and higher order corrections, further 
facilitating the solution process. The major limitation in guidance applications appears to be 
that significant nonlinearities, such as aerodynamic effects must be neglected in the zero 
order problem in order to obtain an analytic solution for the zero order problem, which is 
also nonlinear even in the absence of aerodynamic effects. It turns out in this case that the 
zero order problem is not sufficiently close to the original problem and the solution begins 
to diverge even when a first order correction is attempted (cf. Sec. 3.2b). A second 
drawback which is inherent in any attempt of analysis by model simplification is that a 
significant amount of re-analysis is required when even a minor change in the optimal 
control problem formulation is made. 

Collocation [8, 36] is a general method for obtaining an approximate solution of 
differential equations. It involves choosing simple interpolating functions and enforcing 
the interpolatory constraints at specific points within finite elements to evaluate the 
unknown coefficients. Thus when applied to an optimal control problem, it reduces the 
associated two point boundary value problem to a set of coupled nonlinear algebraic equat- 
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ions. Collocation methods have the advantages that they are simple to use for a wide 
variety of optimization problems, and their accuracy can be improved by increasing the 
number of elements used in the approximation. The major disadvantages are that there is 
no general guarantee that the numerical methods employed will successfully solve the 
nonlinear programming problem undo - all circumstances, and the dimension of the problem 
increases proportionately with the number of elements. 

It is apparent from the above discussion that the advantages of analytical and 
numerical methods are in many respects complementary in the sense that if the advantages 
can be combined in some way, then most of the important disadvantages (from the view- 
point of real time applications) can be reduced. In this chapter, two of possibly many ways 
to obtain such a hybrid methodology are presented, with the potential for use in the 
development of real time optimal guidance algorithms. The first approach uses the method 
of regular expansion to improve upon a collocation solution, thereby reducing the error for 
a given number of elements. The second approach improves upon the first by using both 
regular expansion and analytical methods to identify more intelligent interpolating functions 
in the collocation method, again with the objective of improving the level of accuracy 
without increasing the number of elements. 


4.2 The Method of Collocation 

Collocation is a method for constructing an approximate solution to a set of differ- 
ential equations by using finite elements of polynomials or simple analytic interpolating 
functions. The unknown coefficients are determined by enforcing continuity at the nodes 
and that the time derivatives of the interpolating functions satisfy the differential equations 
at some specified points within each element. We consider an optimization problem with 
unperturbed (ie. g(x, u, t) = 0) dynamics dx/dt = f(x, u, t) and Hamiltonian H = X T f. For 
simplicity, assume a first order polynomial approximation where the derivative constraints 
are enforced at the mid point of each element These constraints can be expressed as: 


Xj-xj— i aH 
tj — tj_2 dX 


t=(tj+tj_ 1 )/2 ; x=(xj+Xj_ 1 )/2 ; X=(A.j+Xj_ 1 )/2 


^j— l _ 3H 
tj-tj-i dx 


t=(t j +tj_ 1 )/2 ; x=(xj+xj_j )/2 ; X=a j +^.j_ 1 )/2 


x(t) = Xj_ 1 +Pj(t-tj_ 1 ) 
X(t) = Xj_ 1 +qj(t-tj_ 1 ) 


44 


; j = 1, ... , N 

» f ® [tj— i , tj] , to = t Q ; t N = t Q + Tq 


(4.1) 


where N is the number of elements. The control is assumed to have been eliminated using 
the optimality condition. In practice, it is more convenient to directly evaluate the nodal 
values (x 0 , X. 0 > •••» x n- ^n) rather than finding the coefficients of the interpolating 
functions. Though higher order polynomials such as Hermite's cubic are generally 
preferred (because of their smoothness properties), we consider a first order representation 
to simplify the presentation, although the approach applies equally well for higher order 
representations. 

4.3 Regular Perturbation Formulation 

A regular perturbation formulation may be introduced by rewriting the actual 
dynamics in the following form: 

x = Pj + e(H x -pj) 
l = qj+e(-H x -qj) 

H u = 0 ; t € [tj_i,tj] (4.2) 

Note that e has again been introduced as a bookkeeping parameter. The justification for 
this step is that if the collocation solution alone accurately approximates the true solution, 
then the second terms in Eq. 4.2 may be regarded as having a small perturbing effects on 
the state and costate derivatives, which is actually zero at the mid points of the elements. If 
the control cannot be eliminated explicitly in the collocation formulation in Eq. 4.1, then an 
analytic portion Il(u, x, X) of the optimality condition (for which it is possible to eliminate 

u) can be extracted such that 

0 = n + e(H u -n) (4.3) 

Note that in the above equations H is the Hamiltonian corresponding to the original system 
without a perturbation parameter. As presented above, a collocation solution may be 
viewed as the zero order solution for the regular perturbation problem formulated in Eq. 
4.2. Also, as will be shown by example in the next section, more intelligent choices of 
interpolating functions can be identified from the necessary conditions, by utilizing to the 
extent possible the analytically tractable portions of the solution. This results in a signi- 
ficant decrease in the computational requirements for a given level of accuracy. 

Now we can apply the perturbation technique described in the Section 3.2 to 
improve the approximate zero order solution from collocation. For the higher order 
problems defined in Eqs. 3.29 - 3.31, we have: 
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A = jgj _ 

1 1 J dx BxBX 

a d 2 H 

A 12 ; = — “ 


t=(tj+i j _ 1 )/2 ; x=(xj+Xj_j)/2 ; X=(Xj+Xj_i>/2 


J BX d Z X 


t=(tj+tj_!)/2 ; x^Xj+Xj.j)^ ;X=(Xj+Xj_i)/2 


A 21: = 


3qj _ a 2 H 


j dx a 2 x 

A dqj b 2 u 

A 22: = — 


t=(tj+tj_j)/2 ;x=(xj+Xj_j)/2 ;X=(Xj+Xj_j)/2 


j dx axdx 

c ij = p j + (t - tj_i)p t j 


t=(tj+tj_!)/2 ; x=(x j +Xj_ 1 )/2 ; X=(A. j +Xj_ 1 )/2 


.Bn 


B 2 H 




C 2 j — qj + Ct — tj_j >q t j 


t=(t j +tj_ 1 )/2 ; x=(x j +Xj_ 1 )/2 ; X=(Xj+X. j _ 1 )/2 


BH 

Bx 


B 2 H 


= (— 


t=(tj+tj_i)/2 ; x=(Xj+Xj_j)/2 ; X=(Xj+Xj_ 1 )/2 


(4.4) 


and for k = 1 , 


Pll- = 


8H 


j BX 


P 21 .=-M 

21 J 8x 


t ; x=x j _ 1 +pj(t-tj_ 1 ) ; X=Xj.,+q j (t-lj. 1 ) P J 
t ; x=Xj_ 1 +p j (t-tj_ 1 ) ; X=X j . 1 +qj(t-tj. 1 )“ q j 


(4.5) 


where all the terms in Eq. 4.4 are constant* within an element, and are evaluated using the 
collocation solution. The matrix Aj is simply the perturbation of the original state and 
costate dynamics evaluated at the constraint point of each element The expression in Eq. 
3.32, which now corresponds to a piecewise constant system matrix Aj, can be written as 


x k(0 

X k (l) 


(t* 


x k (*o) 
^k^o). 



dx + jt n A (t,x) 

L 0 


'Pik' 

. P 2k. 


dx (4. 6) 


* If higher order interpolating polynomials are chosen, the dynamical system will be a 
time varying matrix polynomial with piecewise constant coefficients. 
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and for a constant system matrix Aj, it can be written as 


x k(0 

X k (t) 


= ^A:(^V l) 

J J 


x k(tj-l) + IkQ (t,t j _ 1 )AJ 1 l 
L^k(tj-l)J To A J J J 


Pj 

.‘Ij. 


+A 


-1 


Pi i 

% 


uy 


+ £ ^A:(t^)| 

tj-1 J 


PlkW 

P2k(*> 


dx j t g [ t j , t j 




P ‘j 

L qt jJ 


(4.7) 


where is the state transition matrix and p t j, q^ are defined as in Eq. 4.4. Note that 

is not the same as in Eqs. 3.32 and 3.33 because A is defined differently. The state 
transition matrix here may not have an analytic expression because the zero order solution 
is not necessarily analytic. If this is true, we can solve Eqs. 4.4 and 4.5 using the 
sensitivity functions and superposition property of linear systems. This is done by 
assigning a unit vector for the initial conditions, and numerically integrates the system from 
^ to ^ + T 0 . Thus by changing the position of the non-zero element in the unit vector, the 
sensitivity functions are obtained. This process can be done in parallel for different unit 
vector. 

In the zero order solution, £ in Eq. 4.2 is set to zero, which means that the standard 
collocation constraints in Eq. 4. 1 are employed and an approximate solution is obtained by 
solving the algebraic equations. Then first and higher order corrections may be computed 
by quadrature as explained in the earlier section on regular perturbation. 


4.4 A Duffing's Equation Example 

This investigation is carried out to demonstrate the hybrid approach outlined in the 
preceding section. The example is based on Duffing's equation presented in its first order 
form: 


x = v 


; x(0) = x 0 


v = -x - ax 3 + u ; v(0) = v 0 


(4.8) 


and the objective is to 


mm 

u 


S x x 2 (t f ) + S v v 2 (t f )+ J (1 + — )dt 


(4.9) 


with S x , S v being the weights on the terminal values and tf is free. The problem can be 
converted to the Mayer’s form in Eq. 3.17 (if desired) through the usual method of introdu- 
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cing an additional state equation whose right hand side is the integrand of Eq. 4.9. We 
investigate the problem in different levels of complexity according to how the dynamics of 
the full system are treated, 
a) Level 0 formulation 

This is the degenerate case in which there is an analytic zero order solution, and 
therefore collocation is not required (solely a regular perturbation approach as discussed in 
Sec. 3.2). Let e = a, thus neglecting the hardening effect ax^ in the original problem. The 
necessary conditions are: 


X = V 

O 

X 

II 

© 

X 


3 

v = -x + u - ex 

; v(0) = v 0 


X x = X y + e3A. v x 2 

’ ^x(lf) — 2S x x(tf) 

Xy — X X 

; ^ v (tf) = 2S v v(t f ) 

H u = u + X. v = 0 



|h = A, x v + A, v (-x + u 

“ £X 2 ) + 1 4- U 2 / 2 J 

= 0 
tf 


(4.10) 


The zero order problem (e = 0) is linear and time invariant, and can easily be solved as 


X 

o 

/-"S 

> 

1 


cost 

sint 

(sint- tcost)/2 

-tsint/2 

r x o(t 0 ) i 

vo(0 


-sint 

cost 

tsint/2 

-(sint + tcost)/2 

O 

< 4-* 

o 

> 



0 

0 

cost 

sint 

^xO(^o) 

_A, v g(t)_ 


0 

0 

-sint 

cost 



(4.11) 


where 


t = t - 1„ 


; t 0 ,t e [0,T 0 ] 


(4.12) 


The above state transition matrix is also the state transition matrix for the higher 

order correction. Given the boundary conditions of xo(0) = x 0 , vg(0) = v 0 , ^ x0 (T 0 ) = 
2S x x 0 (T 0 ), X v0 (T 0 ) = 2S v v 0 (T 0 ), the remaining unknowns A, x0 (0), k v0 (0), ^ x o( T o)> 
XvoCTq), Tg can be solved with a Newton's method using Eq. 4. 1 1 and the corresponding 
transversality condition 
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{Hq - *-xO v O ^vO x O ^vO / 2 lj 


= 0 


(4.13) 


From Eq. 3.29, the differential equations governing the higher order correction 
dynamics are 



' x k ’ 


'0 

1 

0 

0 ‘ 

’ x k " 


Vo(0 


>lk(t)* 

d 

v k 


-1 

0 

0 

-1 

v k 

T k 
+ — 

-x 0 (t)-X v0 (t) 

+ 

p 2k(t) 

dt 

^xk 


0 

0 

0 

1 

^xk 

T 0 

A-vO(l) 


P 3k(0 


X v k. 


0 

0 

-1 

0 _ 

_X v k_ 


-A, x o(0 


p 4k(t). 


(4.14) 


with the boundary conditions 

x k (0) = v k (0) = 0 ; ^(Tb) = 2S x x k (T 0 ) ; X vk (T 0 ) = 2S v v k (T 0 ) (4.15) 

In this case, we have for k = 1, 2: 

Pjl = 0 ; P 2 i = — xq ; P 31 = 3X, v qXq , P 44 — 0 

P 12 = w{Ti / T 0 ; P 22 = -(xi + Xyi + x§)T x / T 0 - 3x§ Xl 

P 32 = (^vl + 3^vO x o) T l / T 0 + -^vl x 0 + vO x O x l ’ P 42 = _ ^xl T l I T 0 ( 4 - 16) 


and the transversality conditions: 

{Hj = X x jVo -^vl( x 0 + ^v0) + ^-x0 v l - ^v0( x l +x o)}| 


= 0 


r 2 

|H 2 = X. x 2 v o “ ^v2( x 0 + ^v0) + ^xO v 2 + ^vl x l - 3^ V 0 X 0 X 1 


-A. v i( x i + A. v i + x o) + ^vi /^j 


0 


(4.17) 


which are needed to compute the first and second order corrections by quadrature. The 
results are shown in Figs. 4.1 - 4.4 for S x = S v = 100, and a = 0.4. The first order state 
and costate histories are stored and later retrieved by linear interpolation to construct the 
second order solution. The optimal solution generated using a multiple shooting technique 
is also included for comparison. These results clearly show that the series is not 
convergent, and that the most accurate approximation is obtained using a first order 
solution. If we regard this level of accuracy as insufficient, then the conclusion must be 
that the nonlinear term (ax^) is too large to be neglected in the zero order solution. 
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Lambda- 





b) Level 1 formulation 

This case illustrates the hybrid approach as outlined in the section on collocation, 
using a piecewise linear representation to approximate the states and the costates for the 
zero order solution. The interpolator constraints for an N equally spaced segmentation 
are: 


x 0j ~ X 0j-1 V 0j + V 0j-1 _ 

T 0 / N 2 Pxj 


v 0j V 0j-1 _ x 0j x 0j— 1 ^vOj "*■ ^vOj-1 „/- x 0j x 0j — 1 ^3 _ _ 
rr. /»T ~ ~ I a l Z ) - Pvj 


T 0 /N 

^■xOj ~ ^xOj-1 _ A-yQj + ^vQj-l 


T 0 /N 2 

^vOj ~ ^vOj-1 _ ^xOj ~ ^-xOj-l 

T 0 /N " 2 


l + 3a( X ° i + X ° H ) 2 


4xj 


= q V j 


(4.18) 


with the boundary conditions and transversality condition given by: 

x 00 — x o * v 00 — v o ’ ^xON = ^x x 0N > ^-vON = ^^v v 0N 

a 

^xON v ON + ^vON( _ x ON ~ ^vON ~ ^On) + ^vON / 2 + 1 = 0 (4.19) 


There are 4N+5 equations to solve for the 4N+5 unknowns of xqq, vqq, A^qq, A^qq, ••• * 
X 0N> V 0N’ ^xON’ ^vON> To- Solutions for several values of N are presented in Figs. 4.5 to 
4.8. Note that accuracy improves with increasing N, but at the expense of having to solve 
a large nonlinear system of equations. 

The higher order dynamics in this case are 
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(4.20) 


where 


c = 1 + 3axg 


, b 6a(X v oxo) 


; t = (tj + tj_,)/2 


(4.21) 
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The state transition matrix expression for this case is given in Appendix C. For k - 1,2, 
the forcing function terms in Eq. 4.20 are: 

Pi 1 = v 0 - Pxj ’ p21 = -x 0 — ^v0 " ^0 - Pvj 

P 31 = ^vo(l ’*"3axo) - q x j » P 41 = — ^x 0 “ 9vj 

p 12 = ( v l + v 0 ~ PxjFi / Tq 

P22 = { _cx l “ ^vl - x 0 - ^v0 _ax 0 _ Pvj } T 1 / T 0 - 3ax 0 - X 1 "( 1 + 3ax 0 - c ) x l 

P32 = {cX vl + bxj + X v0 (l + 3axo - q x j)}Ti / T 0 + 6ax 0 -X v l x l + 3a ^v0 - X 1 
+ (1 + 3axo - c)X vl + (6aX v ox 0 - b)xj 

P42 = (— A. x i — X x q — q v j)Ti / T 0 (4.22) 

where 

xo(t) = x 0 j_! + p x j(t - tj_i) ; v 0 (t) = v 0 j_! + p v j(t - tj_i) 

^-xO(t) = ^xOj-l + < lxj(t — *j-l) ; ^vO(0 = ^vOj-1 +t lvj(t — tj-l) (4.23) 

plus the boundary conditions in Eq. 4.14 by replacing x k (0), v k (0), x k (To), v k (To), 
Xx k (T 0 ), ^vkCTo) wi * x k0> v k0’ x kN> v kN> *-xkN< ^vkN- The corresponding expansion of 
the transversality conditions in this case are defined as 

3 3 

0 = A, x in v 0N + ^-xON v IN ~ ^ vlN ( X 0N + ax 0N ) + ^-vON ( — X 1N — ^ vON - ^ON ) 

3 

0 = ^x2N v 0N + ^xON v 2N + ^xlN v lN ~^v2n( x 0N + ax ON) + ^vON( _x 2N 

-X V 2N “ 3axj)N x 2N) + ^vlN( -x lN ~ 3ax 0N x l) - X-vlN / 2 (4.24) 

First and second order corrections are computed for the case where N = 3 is used in 
the zero order collocation solution (note here e is 1.0). The results are shown in Figs. 4.9 - 
4.12. Comparison with the N = 3 results in Figs. 4.5 - 4.8 show that a significant 
improvement in accuracy is achievable without requiring a large number of elements. In 
Figs. 4.9 - 4.12 the second order solution is indistinguishable from the optimal solution. 
The discontinuity in slope (which is a consequence of using first order interpolation 
functions for the collocation solution) is also smoothed as the order of the correction 
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increases. Contrary to Level 0's results, the second order corrections do not diverge due to 
the fact that the nonlinear term has been accounted for in the zero order solution. 



2 3 4 

Time 


Figure 4.5. Level 1 Zero Order Results in x for Different N. 



2 3 4 

Time 


Figure 4.6. Level 1 Zero Order Results in v for Different N. 










Figure 4.1 1. Level 1 Higher Order Results in X* for N=3. 



Figure 4. 12. Level 1 Higher Order Results in Xy for N=3. 


c) Level 2 formulation 

As a second illustration of a hybrid solution approach we retain a portion of the 
dynamics from the necessary conditions to identify a more intelligent interpolating function 
for the hybrid Level 1 formulation. Consider the following simple modification of the 
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regular perturbation formulation for this example: 
x = v 

v = p v j + e|-x - X v - ax 3 - p v j j 
- q x j +e {*-v ( 1 + 3ax ) — Qxj} 

K = ~K (4.25) 

Note that we interpolate only the variables that have nonlinear coupling, and that the result- 
ing interpolation retains more of the dynamics in the original problem than in the Level 1 
formulation. The interpolating functions in this case are: 

x 0(0 = x 0(tj-l) + v 0j— 1 +^Pyj(t-tj-l>]tf- tj_i) 

vo(t) = v 0 j_i+Pvj(t-tj_ 1 ) 

^xO (0 = ^xOj -1 Qxj(t — tj-l) 

^■vo(t) — ^vO(tj-l) — ^xOj -1 fj— l)J(t — fj— l) (4.26) 

Consequently, there are fewer unknowns (2N+5) to be solved and the dynamics retained in 
the formulation should improve the zero order collocation approximation. This allows even 
fewer elements to be used. To evaluate the zero order solution, conditions in Eq. 4.19 are 
enforced by replacing xq^, ^on with xo(tjsj), ^vO^n) from Eq. 4.26, and similarly for 
the first order expressions. The forcing terms for this case are: 

1 — 0 * P 21 — — x 0 — ^vO — ~ Pvj 

P 31 = A, v q( 1 + 3axp) ~ qxj : P 41 = 0 (4.27) 

and the state transition matrix is same as that in Level 1. 

Figures 4.13-4.16 show the zero and first order state solutions for the case N = 2. 
The results show that the zero order solution is dramatically improved especially in the state 
variables in comparison to the zero order solution for N = 3 of the Level 1 formulation in 
Figs. 4.9 and 4. 10. The accuracy of the first order solutions in Figs. 4. 13 to 4. 16 are very 
good and are almost riding on the exact solutions, even though a cruder segmentation has 
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Figure 4. 15. Level 2 Higher Order Results in X x for N=2. 
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Figure 4. 16. Level 2 Higher Order Results in 7^ for N=2. 


d) Level 3 formulation 

In this last demonstration, the Level 2 formulation is further enhanced. All the 
linear terms are retained in the zero order problem, and the nonlinear terms in the v and X y 
dynamics are approximated by piecewise constants. The resultant expressions become: 
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X = V 


v = -x - X, v + p v j + ej-ax 3 - p v j| 
i x =X v + q xj + e{3aX, v x 2 - q xj } 

i v =-K (4-28) 

This is equivalent to the Level 0 problem except for the presence of two additional 
unknown constants. This formulation represents an attempt to make maximum utilization 
of the analytically tractable portion of the solution in selecting the interpolating function for 
the collocation solution in the zero order problem. The zero order solutions in this case are 
also similar to that for the Level 0 case: 

X 0 (t) = (x 0 (tj_l) - Pvj - q x j)cost + v 0 (tj_i)sin t + k x o(tj-l> [sin * - teost] / 2 
-(K vO (tj-1 ) + qxj )t ! sin t / 2 + p v j + q x j 

vq( 0 = -(x 0 (tj_ 1 )-p v j - q x j)sint + v 0 (tj_ 1 )cost + X x0 (tj_ 1 )tsint/2 

-(X, v0 (tj_i) + q x j)[sin t + teost] / 2 

^x()(t) = ^x0(tj-l) cosl + (^vO(Vl> + qxj> sinI 

^ v0 (t) = -^ x0 (t j _i)s i nt + C^vO(tj_i)-t-q x j)cost-q xj (4.29) 


where 


Pvj = -ax 
t = t-tj_l 


; qx J 3a K x2 )>.v((t j +t j _ I )/2);x((t j+ t j _ 1 )/2) 
;t€[t H ,tj] (4.30) 


In this formulation, an efficient way to find the collocation solution is to solve for 
the 2N+5 unknowns of x 0 (0), v 0 (0), X x q( 0), A^O), p vl , q^, p V N, qxN’ tfO usin £ Eq. 
4.29 in Eqs. 4.30 and 4.19. The high order formulations are obtained in the same manner 
as the previous levels and are not repeated here. The zero and first order results using only 
one element are shown in Figs. 4.17 - 4.20. Though the first order results are not as 
accurate as those in Level 2 (because only one element is used), both zero and first order 
solutions are far superior than the Level 0 results (Figs. 4.1- 4.4) which correspond to the 
degenerate case of only one element 
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4.5 Conclusions 

A hybrid analytical/numerical approach for solving optimization problems using 
regular perturbation and collocation methods has been developed. The hybrid approach 
shows that it is possible to significantly improve a collocation solution without increasing 
the number of finite elements. The loss in accuracy that results from using a smaller 
number of finite elements is compensated by the addition of higher order corrections to the 
solution based on regular perturbation theory. Viewed a second way, using collocation to 
solve the zero order problem in a regular perturbation expansion allows more of the 
dynamics to be retained in the zero order solution. It has also shown that further dramatic 
improvements can be achieved by selecting more intelligent interpolating functions which 
are derived from the analytically tractable portions of the necessary conditions. The results 
show important implications in real-time guidance applications which will be demonstrated 
in Chapter 5 on the launch vehicle problem 
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SECTION V 


THE HYBRID APPROACH TO NEAR-OPTIMAL 
LAUNCH VEHICLES GUIDANCE 


This section applies the hybrid analytical/numerical approach of Section 4 to the 
problem defined in Section 2. The feedback guidance approach is based on a piecewise 
nearly analytic zero order solution evaluated using the collocation method. Each piecewise 
representation of the collocation solution obeys a bilinear tangent law for the thrust vector 
angle, which serves as an intelligent interpolating function for the collocation method. The 
zero order solution is then improved through a regular perturbation analysis, wherein the 
neglected dynamics are corrected in the first order term. Wind shear effects and constraints 
are also investigated. 


5.1 Zero Order Solution 

As discussed in Section 4, it is possible to improve a collocation solution by using 
more intelligent interpolating functions than the first order representations in Eq. 4.1. The 
interpolating functions can be derived from analysis of the analytically tractable portions in 
the necessary conditions. In this case if spherical Earth and atmospheric effects are 
neglected then the previous linear tangent law guidance solution results (Sec. 3.2b). 
However, the costate dynamics are poorly represented as either constant or zero. Hence, 
the strategy is to keep the approximation for the state dynamics and use the collocation 
method to improve the representation of the costates (cf. Level 2 and 3 formulation in Sec. 
4.3). This also reduces the number of unknowns by half. Thus instead of using Eq. 4.1 
to interpolate both the states and costates, only the latter are chosen for intetpolation. The 
perturbed collocation formulation in Eq. 4.2 becomes: 


T’(i) _ 50 )a( 1 ) . 

v = vac P ...e-sine- gj } + £ 


m(t) 


(pW - p) A^ sin 0 - sin y + L*^ cos y 


m(t) 


+SP-4+- 

r z r 


T (i) _o( i >A (i) 
u = _vac e_ CQS 0 _ g(i) + £ 
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r = v 


3H 


3H 


^*V *lvj Qvj) > — ^luj tJjjj) » j ~ 1» ••• » N 


.■ . 3H 

^t — Qij ■*" £ ( -v Qrj) 


3v 
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’ 39 


(X v v + A. u u + A, r r) = 0 
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v = Vcosy + W, ;u = Vsiny + W| c 

( fT® -nA® 


H = K 


(T^ -pA^) sin 0-D^ sin y + L® cosy |A e u 2 ^ 

“ * 1 “ 

m(t) 

f (T^ - pA^)cos0 - cosy - siny 


+X. 


r 2 r 


m(t) 


uv 

r 


+ X r v 


(5.1) 


(5.2) 


The terms P^* g§\ g^ are approximations for the average values of the engine nozzle 
back-pressure and the spherical acceleration components for each flight stage. From 
previous investigation it is found that including partial terms for these effects improve the 
approximation, and for the present problem these parameters are chosen as: 

p (1) = p(h 0 )/2 ; gv } =p e /r 2 -u 2 /r 0 igJj^O 

P (2) = 0 ; g$ 2) = gi 1) / 2 ; gjj 2) = 0 (5.3) 

and they are assumed to be updated continuously in closed loop implementation. 

In the following we make use of the analytic portion of the optimality condition in 
Eq. 5.1 to generate the zero order control, by using the form in Eq. 4.3. This amounts to 
regarding the dependence of aerodynamic forces on 0 as a perturbation of the optimality 
condition, which results in the celebrated bilinear tangent law 

<an9o(.) = ^ v0j -‘ +qvi * t ~ 1 i- | | (5.4) 

^■uOj-1 ■** Quj(t tj-l) 


With the above formulation and using the expression in Eq. 5.4 to eliminate the 
control, the zero order solution (e = 0) can be expressed as: 
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vo(t) = v 0 (tj_i) + S]-fc=sinh *(tan((p + rtf) -sinh ^tantp) 

k u> [Vl + A 2 


<P(0 

<P(tj-l) 




(i) 


; t € [tj_j,tj] 


r 

U 0 (t) = u 0 (tj_i) + ^< -j= — =sinh _1 (tan((p + ri))-^sinh _1 (tan(p) 


+ A Z 


9(0 




r 0 (t) = iq( t j_! ) - (A - tan (p)[ f + sinh ^tanty + T^-sinh VtaiKp)] 


sectp-qsinh ^tancp)! 


L VT 

90 j ! 1 ) +[vo(t ^ l, - S ' ) ^ 1(t ' t j- l)_G(t i- l) 


^ vO (0 = ^vO j— 1 + Q vj (* — * j— 1 ) 
^uO(0 = ^uOj— 1 + Quj(t “ tj-l) 
^tO( 0 = ^rOj-1 ^j-l) 


(5.5) 


where 


A — ^ , B — (c v Qvj ^uQuj) / A > C - yCy + Cy B 

D = q v j / A ; E = (c u A - q uj B) / (AC) ; F = T<£ - p (i) A^ 

c v = ^-vOj-l *lvj^j— 1 * c u = ^uOj-1 — *lujO - l ’ c m = +k^tj_j 

_i At + B [ tan -1 (l/ A) , A>0 

(p(t) = tan ( ) ; T| = < , 

C jrc + tan *(1/A) , A <0 


A c m A-f k (i) B . , q uj C 


k (i) C 


c u A q u jB 


c v A - q v jB 

; q = q vj c 


G(t H ) = 


_ ro J q+A j 


sinh 1 (tan((p + ri))-sinh ^tontp) 


k (i) lVT^ 




(5.6) 


The above expressions constitute a set of nonlinear interpolating functions and the zero 
order solution is now expressed in terms of the unknown costate nodal values. To evaluate 
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these values, the collocation constraints on the costate derivatives in Eq. 4. 1 are enforced: 
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l j -t H 
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'dr i 


t— (t j + 1 j_i )/2; X v q — (^ vO j + ^vO j-1 )( 2, 
^•u0 = ^u0j + ^u0j-l) /2; ^r0 = ^r0j + ^r0j-l^ 2 


t— (tj+tj_i)/2; ... ; ^ r o - (^rO j + ^r0 j-1 ^ 2 


t=(tj+tj_!)/2; . . . ; *- r o=(*-rOj + *rOj-l )/2 


(5.7) 


Since more control activity is expected inside the atmosphere, a denser 
segmentation is used for the first stage flight, whereas a 1 -piece segment is sufficient for 
the subsequent more nearly exoatmospheric second stage flight. The total number of 
unknowns to be solved in the zero order problem are 3N+4. Open loop solutions in a 
stationary atmosphere for several increasing values of N are given in Figs. 5.1 to 5.6. The 
segmentation is N-l elements for the first stage flight and one element for the second stage 
flight. Zero order results using only the regular perturbation approach as given in Sec. 3.2 
are also included for comparison. Significant improvements are observed in the costate 
profiles with the hybrid approach because part of the aerodynamic effects are now 
accounted for in the zero order formulation. In particular, note from Figs. 5.4 - 5.6 that the 
zero order solution of Sec. 3.2 amounts to ignoring aerodynamic effects and invoking a 
flat, non rotating Earth approximation. This results in X. u and X r being constant and X v 
being linear in time (see Figs. 5.4 - 5.6), which from Eq. 5.4 gives the linear tangent 
steering law. This largely accounts for the failure of the regular perturbation method when 
aerodynamic effects are included. 
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Figure 5.5. Open Loop Profiles for Various N. 
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Figure 5.6. Open Loop ^ Profiles for Various N. 


5.2 First Order Solution 

In this case, the linear differential equations satisfied by the first order terms have 
the following form: 
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(5.8) 


Complete expressions are given in Appendix D. As explained earlier, the first stage flight 
time is fixed, T = tf - tg, and Tq represents the zero order second stage open flight time. 
Therefore, = 0 for the dynamics describing t < ts, and the second term in Eq. 5.8 is 
dropped for the elements corresponding to this time interval. 

Experience has shown that higher order perturbation corrections are not sensitive to 
using an exact state transition matrix. This behavior is analogous to the practice of using an 
approximate Jacobian to solve nonlinear algebraic equations. So we introduced the 
following approximation to simplify the analysis. The 3x5 lower left comer block of the 
system matrix in Eq. 5.8 represents the effects that second order variations of the 
atmospheric terms have on the costate variables. By neglecting these terms we are able to 
derive an approximate state transition matrix for the first order system: 
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Appendix D details the (0 terms in Eq. 5.9. The lower right hand block in Eq. 5.9 accounts 
for spherical Earth effects on the costate solution, neglected in the zero order solution. As 
will be shown in the numerical results section, this is an important correction for the 
exoatmospheric phase of flight. By successively applying Eq. 4.6 of N times, the 
perturbations at t^ for a first order system with piecewise representation are now given by 


x i(*n)' 
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x l(t 0 )' 
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Jt N-l A N VN Tq 
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P< 2) (X) 

P^X)' 

. J 


j=N-l 

)dx+ £ ^Am^N^N-i) •• 
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ax 


(5.10) 


5.3 Numerical Results 

Figures 5.7 to 5.10 show the closed loop results for the state variables expressed in 
the wind frame coordinates. The control is updated at every second and is held constant 
within each update interval. The total number of elements used in this case is N = 8. Note 
in Fig. 5.10 that jumps in angle of attack occur at about M = 1.3 and M = 2.3. These are 
due to the shadowing effects of the booster which causes the control solution to first follow 
a higher oc profile (to reduce drag) followed by a lower profile to correct the trajectory. 
There is another third small jump at the staging time due to the discontinuous dynamics. 
This figure also shows a major difference between the zero order and first order solution 
for ot during the end of the second stage flight, which is due to the absence of the spherical 
Earth corrections in the zero order solution. Even though a large difference exists between 
the two solutions, the trajectory and the performance index stay very close, and imply that 
the optimal result is insensitive to control variations. 
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Figure 5.9. Closed Loop Altitude Profile for N=8. 



Time (s) 

Figure 5. 10. Closed Loop Angle of Attack Profile for N=8. 

Next, we include the effects of non-stationary atmosphere on the solution. The 
wind profile used is the interpolated mean winter profile for Kennedy Space Center, shown 
in Fig. 2.11, and this profile is accounted for in the guidance solution. From earlier 
investigation it is learned that the performance is not sensitive to control variations, there- 
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Figure 5.11. Convexized First Stage Cp Profile. 
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Figure 5.12. Closed Loop Velocity Profile Under Wind 
and ocq Constraint. 



fore attempt is not made to incorporate the control constraint in the analysis. Instead a hard 
bound on the control is enforced in the simulation, but not in the guidance derivation. The 
bounds in this case are 

-167580 deg Nm' 2 < aq < 167580 deg Nm" 2 (5.1 1) 

They represent the dynamic loading limits on the vehicle. In addition, the first stage Cq 
profile was convexized, as shown in Fig. 5.1 1. This is done to eliminate the objectionable 
jumps in control that are observed in Fig. 5.10, which have negligible effect on the 
performance. The results of the closed loop simulation for the unconstrained case are 
depicted in Figs. 5.12 and 5.15, which show excellent agreement between the first order 
guided solution and the optimal solution. Fig. 5.16 illustrates the effect of the aq 
constraint, which is active only over a minor portion of the trajectory. The performance 
results for this case summarized in Table 5.1. 

Table 5.1. Performance Comparison for ALS Vehicle Guidance. 

optimal lst-order Oth-order 

h(t f ) 148160m 148160.0m 148160.0m 

y(tf) 0° 0.000° -0.001° 

V(t f ) 7858.2ms' 1 7858.20ms' 1 7858.14ms' 1 
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Figure 5.13. Closed Loop Flight-path Angle Profile Under 
Wind and aq Constraint. 
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Figure 5. 16. Closed Loop aq Profile Under Wind 
and aq Constraint. 

In this example, the aerodynamic forces have a major effect in the middle portion of 
the first stage ascent. This is illustrated in Fig. 3.20 which shows the ratios of the lift and 
drag forces to the thrust components. This explains why the regular perturbation analysis 
in Sec. 3.2 is not able to correct for the effect of aerodynamic forces in the first order 
analysis. These forces are simply too large to be treated as perturbation effects, and 
consequently the calculated first order correction diverged. Use of the collocation method 
in forming a zero order solution largely accounts for the aerodynamic effect through the 
mid-element constraints in Eq. 5.7. 

5.4 Remarks on the Numerical Results 

The results show a high level of fidelity and justify the approximation we have 
introduced to obtain the state transition matrix in Eq. 5.9. In particular, the first order 
solution shows significant improvement by correcting for the spherical Earth effects, as 
illustrated in Fig. 5.10. In this case the zero order solution fails to anticipate the sharp 
change in the radial acceleration as the orbital condition is approached, even using a 
continuously updated guess of g v . This results in an excessive pull-up of the vehicle 
during the initial portion of the second stage flight, and is later forced to correct with a large 
negative a to meet the terminal conditions. However, both zero and first order results give 
extremely good orbit injection accuracy without requiring a high rate of control update. 
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Figure 5. 17. A Hypothetical Wind Shear Profile. 



Figure 5. 1 8. Comparison of the Thrust Vector Angle Profiles 
under Wind Shear. 
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The computations for the cases presented here are done on a SPARCstation 1. The 
CPU time needed for a control update ranges from 0.65s for an 8-element case to less than 
0.15s for a 1 -element solution during the second stage flight. The Newton's method with 
Broyden's update of the Jacobian [11] is used in the zero order collocation evaluation and 
the solution converges typically in 4 iterations. It is apparent from the numerical results 
that the first order correction is needed mainly to correct for spherical Earth effects which 
are dominant only in the second stage of flight. Therefore a significant additional savings 
in computation time would have resulted had we computed this correction only for that 
phase. 

5.5 Wind Shear Investigation 

To assess the effectiveness of the hybrid approach against wind shear, we show a 
typical scenario. First, an open loop trajectory using piecewise linear thrust vector angle 
program for the first stage flight, followed by the closed loop hybrid approach guidance for 
the second stage flight is simulated with a hypothetical wind shear (cf. Fig. 5.17). The 
open loop part of the guidance is derived from a linear interpolation of the previous results, 
which is based on the nominal mean wind profile. Second, a guided trajectory using 
closed loop guidance for both the first and second stage flight is simulated. This guided 
solution is assumed to have detected the wind shear, and is therefore included in the 
calculation. To assure structural integrity, both cases are incorporated with the aq 
constraints. The first case represents the approach for present launch vehicle operation, ie. 
an open loop guidance for the endoatmospheric flight using pre-flight atmospheric 
conditions, and compensated by a closed loop guidance for the exoatmospheric flight The 
second case represents the proposed approach for ALS, ie. real-time near optimal guidance. 
Figs. 5.18 - 5.20 compare the ’Open loop' and the ’Guided’ solutions. A point of interest, 
the ’Nominal' solution with the same linear piecewise control program flying under the 
nominal wind condition is also included. The 'Open-loop' solution gives poorer 
performance (cf. Table 5.2). The final time to orbit is 1.13s longer (equivalent to a loss of 
45501bs in payload) than the guided solution, and is also worse than the guided solution 


Table 5.2. Performance Comparison under Wind Shear. 



Guided (1st) 

Guided (0th) 

Open loop 

h(tf) 

148160.0m 

148160.0m 

148160.0m 

7(tf) 

0.000° 

-0.000° 

0.000° 

V(tf) 

7858.20ms* 1 

7858.18ms* 1 

7858.20ms* 1 

-J = tf 

377.287s 

378.243s 

378.413s 
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using only a zero order solution. If the magnitude of the wind shear is further increased by 
22%, then open-loop guidance will result in a catastrophic failure unless the aq limit is 
exceeded. 



0 100 200 300 400 


Time (s) 

Figure 5. 19. Comparison of the otq Profiles under 
Wind Shear. 



Time (s) 

Figure 5.20. Experienced Horizontal Wind Speed for the 
3 Different Simulations. 


82 





SECTION VI 


CONCLUSIONS AND RECOMMENDATIONS 

6.1 Conclusions 

The fundamental problem in treating launch vehicle dynamics by singular 
perturbation methods relates to the inherent large value of longitudinal load factor. As a 
result, the zero order reduced solution gives a very poor approximation. A manifold 
solution was also attempted to account for flight path angle dynamics, but this method also 
fails due to the fact that the dynamics are not separable in the same manner throughout the 
ascent profile. Regular perturbation analysis gives a better solution in the absence of 
aerodynamic forces. However, the approach cannot handle guidance for the atmospheric 
flight phase, which is the main issue of this research. The neglected aerodynamic forces in 
the zero order solution are simply to large to be considered as a perturbation effect. 

A new hybrid approach for the solution of nonlinear problems in optimal control 
has been developed for this application. This approach is hybrid because it combines the 
desirable features of numerical and analytical methods. The numerical method of 
collocation allows a simple formulation for solving a wide variety of optimization 
problems. The disadvantage of requiring a large number of approximation elements and 
solving a large dimension set of algebraic equations are compensated for by the analytical 
approach of regular perturbation. The regular perturbation approach provides higher order 
correction over the collocation solution without increasing the number of approximation 
elements. It can also be used to identify intelligent interpolating functions for the 
collocation solution, which results in a further substantial reduction in the number of finite 
elements needed for a given level of solution accuracy. These attractive features promise an 
enhanced real time capability in the solution of optimal control problems, which has been 
demonstrated in the launch vehicle guidance application. The main results on this problem 
are that a b ilin ear tangent steering law can be employed in all flight phases, including the 
atmospheric phase, and that the collocation solution can be obtained using a small number 
of elements. 

6.2 Recommendations for Future Work 

Many important issues remain for future research, and the following 
recommendations are made in increasing order of complexity: 

Identifying More Intelligent Interpolating Functions - Though the zero order solution is 
capable of handling partial aerodynamic effects, spherical Earth effects were not direcdy 
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incormorpated in the analysis. This leads to a poor representation of the zero order solution 
as the vehicle approaches orbital speed. An investigation should be made to account for 
this effect in the formulation, and a proposed way to set up the problem is to add a constant 
perturbation term, similar to the Level 3 formulation in Sec. 4.3. 

More Acc urate Model - Improvements can be made in the launch vehicle problem by 
considering a more elaborate dynamic model. The rotating Earth effects and a more 
complex propulsion model should be considered. It may be necessary to modify the 
interpolating solution in the collocation approach, depending on the magnitudes of these 
nonlinearities. 

Multi-fligh t Task Requirements - The hybrid guidance approach can be extended to handle 
various flight tasks such as deorbit and rendezvous. These requirements will pose terminal 
constraints on both the downrange and crossrange values, which can be included in a 3-D 
formulation. Such multi-flight task guidance capability would be very useful to manned 
vehicles like the Space Shuttle. 

Constrain ed Problem Analysis - For the launch vehicle problem, it is coincidental that the 
performance is insensitive to control variations, thus allowing the exclusion of the control 
constraint in the analysis. However, it would be useful to complete the hybrid approach to 
include analysis of constrained optimization problems. To address the constrained problem 
requires a guess of the switching structure and a formulation of variable time intervals in 
which the constraint becomes active. 

Launch Vehicles Ra nee Safety Concerns - The range safety issues related to the launch 
vehicle ascent trajectory occur in the form of state constraints. To avoid potential disaster 
or to facilitate the retrieval of reusable boosters, the vehicle may be constrained to fly within 
a narrow corridor of air space. Present methods to handle this type of problem are not 
efficient and rely purely on numerical means. Future study should include a systematic and 
simplified formulation that is tractable by analytic methods such as the hybrid approach. 
Hybrid A pproach with the HJB Expansion - As demonstrated in [14, 37], the regular 
perturbation analysis can be carried out using the Hamilton-Jacobi-Bellman equation. In 
this formulation, the perturbation corrections are not represented by a set of linear O. D. 
E's., and the calculation of the state transition matrix or sensitivity functions are not 
required. Instead the perturbation corrections are evaluated simply by quadrature. 
Proposed future research should include a study of the relationships between these two 
types of formulations (the HJB and the state/costate expansion) and the hybrid approach 
using the HJB expansion, which promises a much simpler and more efficient evaluation of 
the perturbation corrections. 
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Manifold Investigation - The failure of the energy approximation analysis indicates that the 
zero order (e = 0) reduced solution is far from satisfactory and a higher order 
approximation is needed. One distinguishing feature of Manifold Theory is the inclusion of 
e in the manifold condition [24] which considers the fast variable as a function of the 
singular perturbation parameter in addition to the slow variable. Although this approach 
was not successful for this application, due to the varying role of the altitude state, it may 
be highly useful in other nonlinear optimization problems. A drawback in our analysis is 
that we had to numerically experiment to determine a solution close to the manifold. This 
has been accomplished by visual examination of the trajectories in Fig's. 3.5 and 3.6. It 
would be highly desirable to develop an algebraic test for when the initial condition lies on 
the manifold solution. An alternative would be to develop an iterative process that 
converges to the manifold solution. 
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APPENDIX A 
Derivation of Eq. 3.32 


We want to show that 


T 0t. 


'Ci(x)- 

dx = T k 1 to 

*o(0 

Ci(x)_ 

K ryi 

T o 

>o(t). 


(A.l) 


in Eq. 3.32. Let fi - xo(xo,X.q,x); f 2 - X-o(xo,Xo,x), assuming u being eliminated and 
recall that 


Ci=fi+(x-t„)|i 


;C 2 =f 2 +(t-t 0 )^- 

dX 


(A. 2) 


The left hand side of (A. 1) becomes 


ij, 

To, 


f 

d 

f 

V 

\ 

dx 

(*-t 0 ) 

V 

- f 2.. 

/ 


dfj / 3 xq dfj / BXq 


h 

f 2. 


dx (A. 3) 


Using integration by parts on the first term in (A.3), we have 


^j(x-t 0 )ft A (t,x)| 


L f 2J 


A 

A ^ 

~KT-t 0 )(j-fi A (t,x) 
= to T 0 t vdx 


X = 
X 


fl 

f 2j 


dx 


T « t „ 


dfl /dxg dfj / 
df 2 / 3 xq df 2 / d?lo 


fl 

h 


dx 


(A. 4) 


Substituting the state transition matrix property 


d_ 

dx 


^a(L'c) = -^a(Lx) 


dfj / dxo 
df2 / dxo 


dfj / dXg 
df 2 /dX 0 


(A. 5) 


into (A.4), the last two terms cancel and the result is demonstrated. The above state 
transition matrix is also used to derive Eq. 4.7. 
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APPENDIX B 

State Transition Matrix Expression in Eq. 3.44 


where 


The state transition matrix used in the regular perturbation approach in Sec. 3.2 is 



1 

0 

0 

“14 

m (i) 

“l5 

“Sg 1 


0 

1 

0 

“24 

w (») 

“25 

“26 

Q < | ) (t2.tl) = 

l 2 -t l 

0 

1 

“34 

to (i) 

“35 

/v\0) 

“36 

0 

0 

0 

1 

0 

*1 -t 2 
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0 
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0 

0 

0 

0 

1 


°>14 =«14( t 2)-^4 ) (tl) 

©g = 11 1 5 (t 2 ) — *15 (^1 ) 

®16 = *S( l 2) - + t l C0 14 

m (i) - m (i) 

“24 _ ®15 

“25 = *25 to) - *25^1) 

“26 = *26 to) " *S< t l> + 1 1®S 

“34 = 7C §J( t 2)“ 7C 32( t l)“( t 2 -tl)Jti4(tl) 

“35 = *35 to ) ~ *35 (ti) - (t 2 - h )*S (tl ) 

“36 = *36 to) -*36^1)' “to “ t l)*S( t l) +t l“l2 

-sinh -1 [tan(9(t)) - T|)] A sin(6(t)) + cos(6(t)) 
(A 2 +l) 3/2 a 2 +i 


*14 (t) = A i 


*15 (t> = A 


-A sinh 1 [tan(0(t)) - ri)] sin(6 (t)) - Acos(8(t)) | 


(A 2 + l) 3/2 


a 2 + i 
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Q-l 



_ ti f (P + A)sinh 1 [tan(0(t))-ri)] , (p + A)cos(6(t)) + (pA 
*16 w - 7T2"" zm + 


(A 2 *!) 3 ' 2 


a 2 + i 


A = ^vac • B - — 
c u0 k«> ' B -q 

*2 • . - 1 , 


A [ “A 2 sinh ^tan^t))-!])] Asin(0(t)) + cos(6(t)) 

25<t) 1 WPP + a^TI 1 
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(A 2 + l) 3/2 
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A 2 + l I 


(A 2 + l) 3 ' 2 


All the variables are evaluated at the zero order values. 


l)sin(0(t)) 


l)cos(0(t)) 


(B.2) 
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APPENDIX C 

State Transition Matrix Expression of Level 1 Formulation in Sec. 4.3 


The state transition matrix of Level 1 case for the Puffing's example is 




a n 

a 12 

a 13 

a 14 

a 21 

a ll 

— a 14 

a 24 

-ba 24 

b a 23 

a ll 

_a 21 

-ba 2 3 

-bai 3 

-a 12 
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a! x = £cos(Wa) + - — ^-cos(t^) 

a -(3 a-p 

a i 2 = — a ~ C r _ sin(tVa) + — C J^ sin(t^P) 

12 (a-pWa (a-P)# 


a 13 = ' -jg-sina#) 


(a - p)Va 


(a-p)#‘ 

1 
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a 14 = - 2 — cos(tVa) cos(t^/p) 
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a 2l = 
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sinCt^P) 


a 24 = — ^-sin(tVa) + sin(t^/P) 

a-p a-p 


For b < 0: 

6 = (^a/c 2 * -^ -2c)/2 ; 4> = V< c + V"c 2 -b)/2 

2 2 ^ 
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2 2 2 2 
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a 13 

a 21 


~~2 57 sinh (6t) cos(<J>t) + i — y cosh(0t) sin(<bt) 

20(0 2 + <|) 2 ) Y 2<J)(0 2 + <(> 2 ) V 




a 24 = — sinh(0t)cos(«))t)-'^-cosh(0t)sin(«|»t) 

2u 2<J) 


The state transition matrix has the same structure of that in Level 0 for b = 0. 


(C.3) 
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APPENDIX D 

System Matrix and State Transition Matrix Expression 
of the First-orderFormulation in Sec. 5.2 


The terms defined in Eq. 5.8 have the following expressions: 
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The remaining partial derivatives (dq y j/dv, dq v j/du, ... ) are similar to the last three 
expressions in (F.l). 

The approximate state transition matrix in Eq. 5.9 can be obtained by taking the 
partial derivatives of the zero-order solution in Eq. 5.5 with respect to the initial conditions 
( v o(tj-l)» uoOj-lX r o(*j-l)’ ^vo(tj-i), ^uoftj-lX So we have: 
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For example, using the chain rule, C 0 j 4 is given by 

w 3D + ^j^ + dvo^ + 5vo jty + 3^3n_ (D4) 

3D 3c v 3A 3c v 3 £ 3c v 3cp 3c v dr| 3c v 

Symbolic manipulation programs such as Mathematica, MACSYMA can be used to obtain 
the analytic expressions of the above derivatives, and to write the subroutines needed for 
their computation. 
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